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In  this  digital  era  of  our  world,  huge  amounts  of  digital  image  data  are  being  eolleeted  on  a  daily 
basis.  The  eolleeted  image  data  is  being  stored  for  subsequent  proeessing  and  use  in  a  wide  variety 
of  applieations.  For  this  purpose,  it  is  often  important  to  aeeurately  and  preeisely  extraet  relevant 
information  out  of  this  data.  In  eomputer  vision  applieations,  for  instanee,  an  important  goal  is  to 
understand  the  eontents  of  an  image  and  be  able  to  automatieally  gain  an  understanding  of  a  seene, 
implying  an  extraetion  and  reeognition  of  an  objeet.  This  task  is,  however,  greatly  eomplieated  by 
the  aequired  image  data  being  often  noisy,  and  target  objeets  and  baekground  bearing  textural  varia¬ 
tions.  As  a  result,  there  is  a  strong  demand  for  reliable  and  automated  image  proeessing  algorithms, 
for  image  smoothing,  textured  image  segmentation,  objeet  extraetion,  traeking,  and  reeognition. 
The  objeetive  of  this  thesis  is  to  develop  image  proeessing  algorithms  whieh  are  effieient,  statisti- 
eally  robust  and  suffieiently  general,  in  order  to  aeeount  for  noise  and  textural  variations  in  images, 
and  whieh  have  the  ability  to  extraet  and  provide  eompaet  and  useful  deseriptions  of  target  objeets 
in  images,  for  objeet  reeognition  and  traeking  purposes. 

The  main  eontribution  of  the  thesis  is  the  development  of  image  proeessing  algorithms,  whieh 
are  based  on  the  theory  of  eurve  evolution  with  eonneetions  to  information  theory  and  probabil¬ 
ity  theory.  These  eonneetions  form  the  basis  for  extraeting  a  eompaet  objeet  deseription,  in  the 
form  of  a  polygonal  eontour.  One  eontribution  is  the  development  of  a  new  elass  of  eurve  evo¬ 
lution  equations  designed  to  preserve  preseribed  polygonal  struetures  in  an  image  while  removing 
noise.  In  eonjunetion  with  these  flows,  a  loeal  stoehastie  formulation  of  a  well-studied  eurve  evo¬ 
lution  equation,  namely  the  geometric  heat  equation,  provides  an  alternative  mieroseopie  as  well  as 
maeroseopie  view,  whieh  in  turn  led  to  our  proposal  of  vanishing  at  pre-defined  direetions.  Under 
these  flows,  the  limiting  shape  of  a  eurve  is  a  polygon,  pre-speeified  by  the  form  and  the  parameters 
of  the  speeifie  flow.  The  seeond  eontribution  of  the  thesis  is  the  development  of  a  new  aetive  eon- 
tour  model  whieh  merges  the  desirable  polygonal  representation  of  an  objeet  direetly  to  the  image 
segmentation  proeedure  by  adapting  an  information-theoretie  measure  into  an  aetive  eontour  frame¬ 
work  with  an  ultimately  unsupervised  texture  segmentation  goal.  The  polygon-propagating  models 
we  develop  ean  eapture  texture  boundaries  more  reliably  than  the  eontinuous  aetive  eontour  models 
because  the  evolution  of  an  active  polygon  vertex  depends  on  an  overall  speed  function  integrated 
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along  its  two  adjacent  polygon  edges  rather  than  on  pointwise  measurements  along  continuous  con¬ 
tour  points.  In  this  way,  higher-order  statistics  which  provide  more  adapted  information  than  the 
first  and  second-order,  are  captured  through  both  the  nature  of  the  information-theoretic  criterion 
we  utilize,  and  the  nature  of  the  polygon-evolving  ordinary  differential  equations  we  propose.  A 
supplementary  contribution  in  this  sequel  is  a  new  global  polygon  regularizer  algorithm  which  uses 
electrostatics  principles.  The  final  confribufion  of  fhe  fhesis  is  fhe  developmenf  of  a  simple  and 
efficienl  boundary-based  objecf  fracking  algorifhm  well-adapfed  fo  polygonal  objecfs.  This  is  an 
extension  of  fhe  second  confribufion  of  fhe  fhesis,  and  fhe  key  idea  here  is  centered  around  fracking 
a  relafively  few  verfices  fogefher  wifh  fheir  corresponding  edges,  which  in  furn  yields  a  bookkeeping 
simplicity  and  hence  efficiency. 

The  parsimonious  sef  of  feafures  provided  by  fhe  fhree  mefhods  developed  in  Ibis  fhesis  are  use¬ 
ful  for  objecf-based  descripfion  and  recognition  fasks,  and  in  addifion,  may  provide  a  viable  solution 
fo  a  parsimonious,  and  economical  represenfafion  of  large  dafa  sefs  (e.g.  a  confour  represented  by  a 
few  landmarks). 
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Chapter  1 


Introduction 


In  this  digital  era  of  our  world,  huge  amounts  of  digital  image  data  are  being  eolleeted  on  a  daily 
basis.  The  eolleeted  image  data  is  being  stored  for  subsequent  processing  and  use  in  a  wide  variety 
of  applications.  For  surveillance  applications  for  instance,  the  amount  of  remotely  sensed  imagery 
collected  daily  by  space-borne  or  airborne  systems  is  on  the  order  of  terra  bytes.  It  is  often  important 
to  accurately  and  precisely  extract  relevant  information  out  of  this  data  for  a  variety  of  applications. 
It  is,  however,  impossible  to  manually  carry  out  this  processing  and  analysis  wholly  by  human  op¬ 
erators.  The  acquired  image  data  is  often  noisy,  and  target  objects  and  background  bear  significant 
textural  variations.  It  may  also  be  desirable  to  track  features  or  objects  in  an  image  through  time  to 
obtain  a  dynamic  analysis  of  the  scene.  In  computer  vision  applications,  for  instance,  an  important 
goal  is  to  understand  the  contents  of  an  image  and  be  able  to  automatically  gain  an  understanding 
of  the  scene  and  the  surroundings,  implying  an  extraction  and  recognition  of  an  object.  As  a  re¬ 
sult,  there  is  a  strong  demand  for  reliable  and  automated  image  processing  algorithms,  for  image 
smoothing,  textured  image  segmentation,  object  tracking  in  video  sequences,  and  object  extraction 
and  recognition. 

The  broad  objective  of  this  thesis  is  that  of  developing  image  processing  algorithms  which  are 
efficient,  in  the  sense  of  ease  of  computations,  fast,  statistically  robust,  in  the  sense  of  being  resilient 
to  noise,  statistically  significant  and  meaningful,  in  the  sense  of  accounting  for  textural  variations  in 
images,  with  an  ability  to  extract  and  provide  compact  and  useful  descriptions  of  target  objects  in 
images,  for  object  recognition  and  tracking  purposes. 

We  next  give  the  motivations  for  the  image  processing  algorithms  we  developed,  which  con¬ 
stitute  the  core  contributions  of  this  thesis,  and  then  describe  the  organization  and  summary  of  the 
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thesis. 


1.1  Thesis  Motivations  and  Contributions 

Curve  evolution  teehniques  for  image  proeessing  involve  a  propagation  or  deformation  of  a  eurve  via 
certain  partial  differential  equations  (PDEs).  Curve  evolution  techniques  are  applied  to  a  variety  of 
problems,  such  as  image  filtering,  smoothing,  image  segmentation,  object  tracking,  shape  analysis, 
and  morphological  operations.  The  main  contributions  of  this  thesis  are,  first  the  development  of 
a  new  class  of  curve  evolutions  as  nonlinear  filters  for  curve  and  image  smoothing,  second  the 
development  of  a  new  class  of  curve  evolutions,  where  the  curve  takes  the  form  of  a  polygon  for 
image  segmentation,  and  which  makes  use  of  an  information-theoretic  measure  adapted  to  texture 
segmentation.  A  third  contribution  of  the  thesis  is  the  extension  of  the  second  algorithm  that  is 
developed  in  Chapter  4  to  object  tracking  in  video  sequences.  Potential  applications  of  the  proposed 
algorithms  in  shape  or  object  recognition  are  also  alluded  to. 

The  motivation  and  contribution  of  each  of  these  algorithms  are  described  in  the  next  three 
subsections. 

1.1.1  A  New  Class  of  Curve  Evolutions  for  Nonlinear  Filtering 

Curve  and  image  evolution  techniques  have  emerged  in  recent  years  as  important  applications  of 
PDEs  to  nonlinear  curve  and  image  filtering.  The  well-known  low-pass  filtering  of  an  image  by  a 
Gaussian  function  has  been  shown  by  Witkin  [147]  to  be  equivalent  to  evolving  an  image  with  a 
heat  equation,  which  is  a  linear  PDE  also  known  as  a  linear  diffusion.  This  in  turn  lead  to  vari¬ 
ous  developments  in  nonlinear  filtering  of  images  through  numerous  PDEs  and  nonlinear  diffusion 
techniques  [6, 110, 117, 122, 136].  In  addition  to  these  developments  in  image  filtering,  the  advance 
of  curve  evolution  techniques  may  be  considered  as  an  implication  of  the  desire  to  develop  shape- 
based  techniques,  which  is  closer  in  spirit  to  providing  object  level  knowledge  in  image  analysis. 
Image  evolution  equations  operate  on  a  pixel-based  knowledge,  whereas  curve  evolution  equations 
operate  on  individual  level  curves  of  an  image,  hence  at  a  higher  level  than  the  former.  Shape  infor¬ 
mation  of  an  object  in  an  image,  may  be  reminiscent  of  a  curve  extracted  from  an  image,  thus  for 
an  object-oriented  filtering,  curve  evolutions  may  be  more  appropriate. 

Much  of  the  research  in  curve  evolution  theory  has  centered  around  the  so-called  geometric  heat 
equation,  which  is  well-known  for  its  smoothing  properties.  Indeed,  any  curve  evolved  through  the 
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geometric  heat  equation  is  circularized  [54].  One  way  of  overcoming  this  effect  is  by  using  prior 
information  on  the  preferred  directions  of  important  features  of  objects  in  an  image.  Polygonal 
structures  are  ubiquitous  in  images  of  man-made  objects  such  as  buildings,  roads  and  so  on,  which 
contain  many  straight  lines,  often  oriented  in  particular  directions. 

The  first  contribution  of  this  thesis  is  the  development  of  curve  evolution  techniques  which 
constitute  a  new  class  of  nonlinear  filters  aimed  for  smoothing  along  salient  lines  of  an  image.  These 
filters  which  act  on  level  curves  preserve  sharp  corners,  and  hence  preserve  the  polygonal  structures 
present  in  an  image.  This  is  achieved  by  a  modification,  indeed  a  directional  generalization,  of 
the  celebrated  geometric  heat  equation.  The  designed  equations  or  flows  preserve  certain  classes 
of  features  in  a  curve.  This,  followed  a  local  stochastic  formulation  of  the  geometric  heat  flow, 
leading  to  a  new  macroscopic  view  of  this  equation,  which  is  later  further  specified  to  vanish  at 
pre-defined  directions.  The  limiting  shape  resulting  from  each  flow  that  belongs  to  the  new  class 
of  curve  evolving  flows  is  a  polygon,  pre-specified  by  the  form  and  the  parameters  of  the  specific 
flow.  When  applied  to  an  image,  a  selective  nonlinear  filtering  along  the  salient  lines  in  the  image 
is  achieved  slowing  down  the  effects  of  geometric  diffusion  across  important  structural  directions. 
This  new  class  of  filters  are  presented  in  Chapter  3. 

We  have  also  developed  a  variant  of  this  contribution,  for  smoothing  a  group  of  level  curves 
of  synthetic  aperture  radar  (SAR)  imagery,  which  is  particularly  robust  to  speckle  noise  typically 
present  in  SAR  images.  These  curve  evolutions,  which  lead  to  a  simple  segmentation  of  SAR 
images  and  target  recognition  on  the  extracted  silhouettes,  are  presented  in  [140, 142]. 

1.1.2  A  Polygon  Evolution  Approach  to  Image  and  Texture  Segmentation 

Image  segmentation  has  received  a  lot  of  attention  through  the  years,  and  a  vast  array  of  differ¬ 
ent  approaches  have  been  proposed,  such  as  thresholding  and  local  filtering  [22,97, 112, 134],  re¬ 
gion  growing  [100],  snakes,  and  active  contours  [25,34,70,72,93],  and  global  energy  minimiz¬ 
ing  techniques  using  various  perspectives  such  as  Bayesian,  and  minimum  description  length  [17, 
51,87, 102].  Following  the  pioneering  snake  methodology  [70],  curve  evolution  approaches,  so- 
called  active  contours,  have  been  popular  in  image  segmentation.  The  key  idea  in  the  active  con¬ 
tour  framework  is  the  construction  of  an  energy  functional  for  the  active  contour,  and  its  min¬ 
imization  through  the  gradient  descent  equations  that  propagate  the  active  contour.  Two  main 
categories  of  active  contour  methods  are  given  by;  geometric  (also  called  geodesic)  active  con¬ 
tours  [23,25,72,93],  where  the  curve  evolution  is  based  on  edge  information,  and  region-based 
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active  contours  [29, 115, 121, 150],  where  the  curve  evolution  is  based  on  the  region  information 
inside  and  outside  the  curve. 

The  region-based  active  contour  models  proved  to  be  more  robust  to  noise  conditions  when 
compared  to  edge-based  active  contour  models.  The  region-based  model  usually  utilizes  simple 
region  statistics  such  as  means  and  variances,  hence  can  not  account  for  higher  order  nature  of 
the  textural  characteristics  of  image  regions.  In  addition,  the  object  delineation  by  active  contour 
methods,  results  in  a  contour  representation  which  still  requires  a  substantial  amount  of  data  to  be 
stored  for  subsequent  multimedia  applications  such  as  visual  information  retrieval  from  databases. 
Polygonal  approximations  of  the  extracted  continuous  curves  are  required  to  reduce  the  amount  of 
data  since  polygons  are  powerful  approximators  of  shapes  for  use  in  later  recognition  stages  such 
as  shape  matching  and  shape  coding. 

The  second  contribution  of  this  thesis  is  the  development  of  a  new  active  contour  model  which 
nicely  ties  the  desirable  polygonal  representation  of  an  object  directly  to  the  image  segmentation 
process  by  including  an  information-theoretic  measure  into  the  active  contour  framework  with  an 
unsupervised  texture  segmentation  goal.  The  polygon-propagating  models  we  develop  can  robustly 
capture  texture  boundaries  relative  to  the  continuous  active  contour  models  as  the  evolution  of  an 
active  polygon  vertex  depends  on  an  overall  speed  function  integrated  along  its  two  adjacent  poly¬ 
gon  edges  rather  than  pointwise  measurements  along  continuous  contour  points.  In  this  way,  more 
higher-order  statistics  than  the  first  and  second-order  are  rationally  captured  through  the  proposed 
information-theoretic  measure  we  utilize,  and  the  nature  of  the  polygon-evolving  ordinary  differen¬ 
tial  equations  we  derive.  This  new  variational  texture  segmentation  model,  is  unsupervised  since 
no  prior  knowledge  on  the  textural  properties  of  image  regions  is  used,  and  will  be  described  in 
Chapter  4. 

A  by-product,  nevertheless  necessary,  contribution  in  this  sequel  is  a  new  polygon  regularizer 
algorithm  which  uses  electrostatics  principles.  This  is  a  global  regularizer,  and  is  more  consistent 
in  preserving  local  features  such  as  corners  than  a  local  polygon  regularization,  as  is  explained  in 
Section  4.5. 

1.1.3  A  Polygon  Propagation  Approach  to  Video  Object  Tracking 

Video  sequences,  i.e.  time-varying  image  sequences,  provide  additional  information  on  how  scenes 
and  objects  change  over  time  when  compared  to  still  images.  The  problem  of  tracking  moving 
objects  has  a  variety  of  applications  in  video  surveillance,  traffic  moniforing,  video  coding,  and 
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robotics.  Object  tracking  methods,  may  be  classified  into  two  categories  according  to  the  type  of 
information  they  use.  Boundary-based  methods,  which  utilize  the  boundary  information  along  the 
object’s  contour,  make  use  of  snake  models  [12,  67, 70,  89],  or  geometric  active  contour  models  [24, 
114],  and  usually  constrain  motion  by  certain  motion  models  such  as  rigid  or  affine.  Region-based 
methods  [9, 10, 145]  segment  an  image  sequence  into  regions  (with  different  motions),  which  are 
matched  to  estimate  motion.  The  cost  of  matching  regions  significantly  increases  the  computational 
burden  of  these  techniques. 

The  third  contribution  of  this  thesis  is  the  development  of  a  simple  and  efficient  boundary-based 
tracking  algorithm  well-adapted  to  polygonal  objects.  We  build  on  the  insight  gained  from  the 
the  second  contribution  of  the  thesis,  namely  the  active  polygon  framework,  to  extend  it  to  track 
polygon  vertices  in  time-varying  images.  The  key  idea  here  is  centered  around  tracking  a  relatively 
few  vertices  together  with  their  corresponding  edges,  which  in  turn  yield  a  simplicity  and  efficiency 
in  bookkeeping.  This  objecf  tracking  method,  together  with  an  experimental  study  of  applying 
active  polygons  to  an  object  recognition  scenario,  are  presented  in  Chapter  5. 

1.1.4  Connections  Among  the  Contributions 

The  three  main  contributions  of  this  thesis  may  be  cast  within  a  unified  objective  of  extracting  a 
compact  object  description,  in  the  form  of  a  polygonal  contour,  which  leads  to  an  efficient  repre¬ 
sentation  of  an  object  crucial  to  subsequent  computer  vision  applications.  The  first  contribution,  in 
Chapter  3,  aims  at  removing  unwanted  perturbations  on  curves  while  preserving  salient  features.  It 
drives  a  curve  (a  level  curve  of  an  image),  which  is  assumed  to  contain  shape  information  of  an  ob¬ 
ject  in  an  image,  towards  a  polygon  by  straightening  the  curve  out.  Inspired  by  these  filters  designed 
for  image  smoothing,  we  proceed  to  a  direct  polygonal  description  of  an  object  in  an  image  with  a 
resulting  segmentation,  which  is  the  second  contribution  of  this  thesis,  in  Chapter  4.  This  compact 
description  of  an  extracted  object  by  a  handful  of  vertices,  in  turn  leads  to  the  idea  of  tracking  these 
features  in  a  time- varying  image  sequence,  as  elaborated  in  Chapter  5.  This  contribution  may  also 
be  viewed  as  a  generalization  or  extension  of  that  of  Chapter  4. 

The  parsimonious  set  of  features  provided  by  the  three  methods  developed  in  this  thesis,  are  use¬ 
ful  for  object-based  description  and  recognition  tasks,  and  in  addition,  may  provide  a  viable  solution 
to  a  parsimonious,  and  economical  representation  of  large  data  sets  (e.g.  a  contour  represented  by  a 
few  points). 
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1.2  Thesis  Summary  and  Organization 


In  Chapter  2  we  provide  a  baekground  on  the  teehniques  and  eoneepts  that  are  relevant  throughout 
this  thesis.  A  review  on  image  and  eurve  evolution  teehniques  is  presented  in  Seetion  2.2,  and  the 
eorresponding  numerieal  implementation  teehnique  is  given  in  Seetion  2.2.3.  A  brief  review  of  the 
literature  on  segmentation  methods  is  presented  in  Seetion  2.3.  Snakes,  and  aetive  eontour  models 
that  are  related  to  the  algorithms  developed  in  this  thesis  are  introdueed  in  Seetion  2.3.3. 

Chapter  3  deseribes  a  new  elass  of  eurve  evolutions  developed  for  feature-preserving  eurve 
and  image  filtering  with  a  prior  knowledge  on  the  salient  line  direetions  of  objeets  in  an  image. 
The  theory  of  stoehastic  differential  equations  (SDEs)  is  briefly  introdueed  in  Seetion  3.3.1,  and  a 
formulation  of  the  eelebrated  geometrie  heat  equation  by  a  loeal  SDE,  and  thus  obtaining  a  new 
mieroseopie/maeroseopie  view  of  this  equation,  is  derived  in  Seetion  3.3.2.  The  insight  led  to  de¬ 
signing  a  new  elass  of  nonlinear  filters  in  the  form  of  diffusions  that  vanish  at  pre-defined  direetions 
as  explained  in  Seetion  3.4.  These  diffusions  on  eurves  yield  a  limiting  polygon  shape  preseribed  by 
the  parameters  of  the  flow.  We  have  applied  this  algorithm  to  smoothing  of  struetures  along  known 
orientation  of  salient  lines  in  an  image  while  preserving  important  features.  We  suggest  further  ap- 
plieations  of  the  proposed  flows  in  Seetion  3.5  sueh  as  a  shape  morphing  applieation  in  eomputer 
graphies  and  a  shape  reeognition  seenario. 

In  Chapter  4,  we  present  a  new  elass  of  variational  aetive  eontour  models  for  an  unsupervised 
texture  segmentation  problem.  A  brief  literature  review  on  texture  analysis  and  segmentation  is 
given  in  Seetion  4.2.1.  Through  a  eombination  of  a  novel  polygon  evolution  model  in  Seetion  4.3, 
and  an  information-theoretie  eriterion  adapted  as  an  energy  funetional  for  the  aetive  polygons  in 
Seetion  4.4,  a  robust  texture  segmentation  algorithm  is  developed  as  validated  through  extensive 
simulation  results  provided  at  the  end  of  Chapter  4.  A  generalization  of  the  proposed  aetive  eontour 
model  to  evolution  of  multiple  aetive  eontours  is  also  given.  Seetion  4.5  presents  a  novel  global 
polygon  regularizer  idea  with  a  goal  of  avoiding  degeneraey  during  propagation  of  the  aetive  poly¬ 
gons. 

Chapter  5  extends  the  aetive  polygons  aeting  on  a  single  image  to  time-varying  image  sequenees, 
and  presents  a  video  objeet  traeking  method  in  Seetion  5.1,  with  an  applieation  to  traeking  targets  in 
infra  red  (IR)  image  sequenees.  A  brief  overview  on  objeet  traeking  methods,  and  motion  estimation 
methods  are  also  given  respeetively  in  Seetion  5.1.1,  and  Seetion  5.1.2.  We  further  demonstrate  the 
utility  of  the  aetive  polygons  with  an  experimental  study  on  an  objeet  reeognition  applieation  in 
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Section  5.2.  This  may  be  considered  as  part  of  future  research  based  on  the  framework  set  in  this 
thesis. 

Conclusions  and  possible  directions  of  future  research  are  presented  in  Chapter  6. 
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Chapter  2 


Preliminaries 


This  chapter  focuses  on  the  background  of  techniques  and  concepts  that  are  of  relevance  throughout 
this  thesis. 


2.1  Notation 

A  digital  image  to  be  processed  is  a  2-Dimensional  (2-D)  function  denoted  by  7,  7  :  — )■  M,  where 

C  is  the  domain  of  the  function.  Processing  a  function  7°(rE,  y),  which  depends  on  two  spatial 
variables,  rr  G  M,  and  y  G  M,  via  a  partial  differential  equation  (PDE)  takes  the  form; 


It 

1(0, x,y) 


•1^(1,  Ix,  ly,  Ixx,  Ixy,  lyy) 

I%x,y). 


(2.1) 


Here  t  is  called  the  time  or  scale.  If  denotes  the  partial  derivative  with  respect  to  (w.r.t.)  t  (sometimes 
shown  as  ^),  and  partial  derivatives  w.r.t.  spatial  variables  are  shown  inside  the  operator  A  on  the 
right  hand  side.  The  solution  7(f,  x,  y)  is  referred  to  as  a  scale  space  for  0  <  f  <  oo.  If  ,4  is  a  linear 
(nonlinear)  operator,  the  scale  space  is  called  linear  (nonlinear). 

Definitions  of  some  operators  that  are  commonly  invoked  in  PDE’s  often  used  in  computer 
vision  are  given  next.  The  gradient  of  7 (f,  x,  y)  is  a  two-dimensional  vector  (2-D  vector)  defined  as 


where  the  superseript  T  denotes  the  transpose  of  a  veetor.  The  L2  norm  of  the  gradient  is  given  by 


llv7||tyri  +  /2. 


(23) 


Veetor  fields  arise  when  the  gradient  operator  V  is  applied  to  a  sealar  funetion  sueh  as  I{x,  y).  The 
divergence  of  any  veetor  field  {I{x,  y),  J{x,  y))^  is 


J 


^  I  +  J 

—  ■‘■X  ' 


(2.4) 


The  Laplacian  operator  aefing  on  a  2-D  funelion  is  defined  as  A7  =  Ixx  +  lyy  It  is  elear  that 
the  Laplaeian  ean  also  be  written  as  A7  =  V  •  V7.  Throughout  the  thesis,  an  inner  produet  (dot 
produet)  is  denoted  by  either  <  •  >,  or  by  a  dot  •  .  The  area  of  a  region  R  is  denoted  by  |7i!|.  A 

boldfaee  letter  is  used  to  denote  both  a  veetor  and  a  matrix. 

Some  baekground  on  eurve  and  image  evolution  teehniques,  and  an  overview  of  image  segmen¬ 
tation  teehniques  are  presented  in  the  following  seetions. 


2.2  Image  and  Curve  Evolution  Techniques 

Obtaining  a  family  of  images  (eurves)  from  an  initial  image  (eurve)  through  a  PDE  is  referred  to  as 
an  evolution  of  the  image  (eurve)  through  time  t.  It  is  equivalent  to  the  seale  spaee  eoneept  defined 
in  fhe  previous  seefion.  We  look  af  some  well-known  evolufion  equations  for  images  and  eurves 
nexf. 

2.2.1  Image  Evolutions 

If  is  known  fhaf  a  low  pass  Gaussian  filler  from  signal  proeessing  may  be  implemented  by  evolving 
fhe  intensifies  of  an  image  R{x,  y)  via  fhe  linear  heal  equation  [147], 

I{^:X,y)  =  I%x,y), 

It{t,x,y)  =  V  •  (V7(f,a;,y)),  f  >  0.  (2.5) 

The  solution  to  Ihis  equation  yields  a  parameterized  family  of  new  images  I{t,x,y),  where  fhe 
image  al  eaeh  time  f  >  0  is  equivalent  to  the  original  image  Io{x,  y)  =  7(0,  x,  y)  eonvolved  with  a 
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Gaussian  filter  of  variance  2t.  This  equivalence  gives  rise  to  a  natural  generalization  of  the  low  pass 
filter  using  nonlinear  diffusion. 

Nonlinear  diffusion  has  a  distinct  advantage  in  image  processing  over  linear  diffusion  in  that 
it  may  be  allowed  to  handle  anisotropies  (giving  rise  to  the  so-called  anisotropic  dijfusion)  in  an 
image. 

A  popular  approach  to  anisotropic  diffusion  is  based  upon  models  first  introduced  by  Perona  and 
Malik  in  [117].  Since  then,  these  models  have  received  a  tremendous  amount  of  attention,  as  have 
the  models  based  upon  curve  evolution  theory.  Perona  and  Malik  extended  the  linear  heat  equation 
by  considering  diffusion  coefficients  which  vary  with  the  strength  of  the  gradient  at  different  points 
of  an  image.  This  leads  to  PDE’s  of  the  form  If  =  V  ■  (p(||V/||)V/),  where  p  :  M  — )■  M  is  typically 
a  monotonically  decreasing  function  which  suppresses  diffusion  where  the  gradient  is  high  (near  an 
edge). 

Nonlinear  diffusion  is  particularly  important  when  salient  image  features  are  of  interest.  For 
example,  when  the  preservation  of  sharp  edges  is  important,  it  is  natural  to  consider  an  anisotropic 
model  which  diffuses  an  image  only  along  the  local  direction  of  its  edges.  One  such  approach  is  to 
consider  an  image  I{x,  y)  as  a  collection  of  iso-intensity  contours,  or  level  curves,  and  to  note  that 
at  an  edge  point,  the  direction  of  the  edge  corresponds  to  the  tangent  of  the  iso-intensity  contour 
running  through  that  point.  Let  y  denote  the  direction  normal  to  the  level  curve  through  a  given 
point  (the  gradient  direction),  and  let  ^  denote  the  tangent  direction. 

We  may  write  these  directions  in  terms  of  the  first  derivatives  of  the  image  as 

^  _  {Ixi  ly)  ^  _  {~Iyi  ^x) 

Since  these  constitute  orthogonal  directions,  we  may  exploit  the  rotational  invariance  of  the  Lapla- 
cian  operator  and  re-write  the  linear  heat  equation  in  terms  of  these  two  variables: 

=  V  •  (V/)  =  Irjn 

where  Inn  and  %  denote  the  second-order  directional  derivatives  in  the  directions  of  rj  and  ^  re- 
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spectively.  One  may  then  derive  the  following  expressions 


Iriri  — 


%  = 


pT  _i_2rrr 


^y-^xy 

i^  +  il 


/■2/-  —2TTT  -\- 

■‘■y’-xx  ^-^xJ-yJ-xy  ^  J-x-^yy 


(2.6) 

(2.7) 


By  subtracting  the  normal  diffusion  component  (2.6)  from  the  linear  heat  equation,  which  diffuses 
isotropically,  we  obtain  the  following  anisotropic  model,  which  diffuses  along  the  boundaries  of 
image  features  but  not  across  them 


It  =  he  = 


/■2r  —  2T  T  T  T‘^T 

■‘■yi-xx  ^J-xJ-yJ-xy  ^  ^x-‘-yy 


(2.8) 


We  may  obtain  this  same  equation  in  a  completely  different  and  much  more  geometric  manner  by 
specifying  the  evolution  of  each  level  curve  in  the  image  as  seen  in  the  next  section. 


2.2.2  Curve  Evolutions 

Let  us  denote  a  family  of  smooth  curves  C  {p,  t)  =  {X{p,  t),  yip,  t)),  which  is  a  mapping  from 
X  C  M  X  [0,  T]  — )■  where  p  G  X  is  a  parameter  along  the  curve,  and  t  parameterizes  the  family 
of  curves.  Denote  the  tangent  vector  to  the  curve  at  p  by  T  =  C  =  and  the  normal 

vector  to  the  curve  at  p  by  AT  =  [—y^  X').  We  consider  regular  curves  whose  tangent  is  never 
zero  (C'(p)  /  0  for  allp  G  X). 

Given  p  G  X,  the  arclength  of  a  regular  parameterized  curve  C  from  the  point  Po  is  by  definition 

s{p)=  r \\C'{p)\\dp,  (2.9) 

Jpo 

where  ||C'(p)||  =  ^JiX'{p)Y  +  {y'ipy  is  the  length  of  the  vector  C  [41].  Hence  ds/dp  = 

\\C'ip)\\. 

The  most  general  deformation  of  a  planar  curve  C  °  is  given  by 


dC 

~dt 

C(p,0) 


a{p,t)T  +Pip,t)N , 

Chp). 


(2.10) 

(2.11) 


It  can  be  shown  that  by  a  reparameterization  of  points,  the  first  equation  above  can  be  reduced  to  [75] 
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=  P{p,t)N  .  Considering  curve  evolutions  which  depend  only  on  the  curvature  function  k 
of  the  curve,  f3  can  be  written  as  F{k),  and  a  local  deformation  as  a  function  of  curvature  may  be 
written  as 

Ct  =  F{K)N.  (2.12) 

The  curvature  function  is  the  second  derivative  of  the  curve  C ,  in  the  direction  of  the  unit  normal 
N .  If  the  curve  is  parameterized  by  its  arclength  parameter  s,  the  second  derivative  is  given  by 
C  ss  =  .  The  following  flow 

Ct  =  Css  =  nN,  (2.13) 

referred  to  as  the  Geometric  Heat  Equation  (GHE),  is  well  known  for  its  smoothing  properties.  It 
has  been  shown  by  Grayson  [54]  that  any  closed,  embedded  curve  evolving  according  to  (2.13)  will 
convexify  and  smoothly  shrink  to  a  single  point  in  finite  time,  becoming  more  and  more  circular 
along  the  way.  This  flow  is  also  referred  to  as  the  curve  shortening  flow  since  it  corresponds  to  the 
gradient  (descent)  evolution  of  the  arclength  functional.  See  [74-76]  for  a  more  extensive  discussion 
of  the  many  properties  associated  with  this  flow.  Since  the  evolution  speed  is  a  function  of  the 
curvature  at  each  point  on  a  curve,  this  flow  gives  rise  to  a  Euclidean  invariant  scale  space  (see  [5, 
6, 147])  in  which  finer  features  are  removed  first,  followed  by  coarser  features,  as  the  curve  evolves. 
A  related  flow,  based  upon  the  affine  geometry  of  the  curve,  is  given  by  and  shares 

many  of  the  same  properties  as  the  curve  shortening  flow  but  gives  rise  to  a  more  general  affine 
invariant  scale  space  (see  [5, 124, 125]). 

2.2.3  Level  Set  Method 

A  new  and  efficient  method  for  evolving  a  single  iso-intensity  contour  has  recently  been  proposed, 
referred  to  as  a  level  set  method  [126].  The  parameterized  curve  C  {p,  t)  is  embedded  into  a  surface, 
which  is  called  the  level  set  function  $(rE,  y,  f)  :  x  [0,  T]  i-)-  M.  The  curve  C  is  the  zero-level 

set  of  this  function  ^{x,y,t): 


C  =  {{x,y)  :  ^{x,y,t)  =0}.  (2.14) 

The  evolution  equation  for  ^{x,y,t)  is  derived  from  the  constraint  that  at  any  time  t,  we  should 
have  $(C  {t),t)  =  ^{X{t),y{t),t)  =  0,  and  differentiating  this  constraint  with  respect  to  t,  we 
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get 


+  ^x^t  +  ^yyt 

=  0, 

(2.15) 

%)  ■  ytf 

=  0, 

(2.16) 

-h  V  $  •  C  t 

=  0. 

(2.17) 

Substituting  the  general  form  of  the  curve  evolution  equation  Eq.(2.12),  which  depends  on  local 
geometry  of  the  curve,  into  Eq.(2.17)  above  yields, 


+  V$  •  F(K)Ar  =  0.  (2.18) 

Noting  that  the  outward  unit  normal  vector  can  be  written  as,  N  =  V$/||V$||,  an  evolution 
equation  for  $  is  given  by 

= -F(k)||V$||.  (2.19) 

Thus,  the  curve  C  evolving  according  to  Eqn  (2.12)  can  be  obtained  by  the  zero-level  set  of  the 
function  $  which  evolves  according  to  Eqn  (2.19).  The  selection  of  the  speed  function  F{k)  has 
been  a  subject  of  research  [126].  The  simplest  form  where  F{k)  =  —k  results  in 


=  K 


|V$| 


(2.20) 


Note  that  the  unit  normal  and  the  curvature  of  a  level  curve  may  be  expressed  as  iV  =  |j^§||  and 
K  =  V  •  ^  ^ .  This  allows  us  to  rewrite  the  above  equation  completely  in  terms  of  $  and  its 

derivatives. 


\  xy  +  yy 


(2.21) 


giving  us  a  PDE  which  is  identical  to  (2.8),  and  is  also  referred  to  as  the  geometric  heat  equation 
since  it  is  a  result  of  applying  the  previous  geometric  heat  equation  (2.13)  to  the  zero-level  curve  of 
the  level  set  function  $. 

Eor  curve  propagation,  the  level  set  methods  constitute  an  Eulerian  approach  in  which  the  un¬ 
derlying  coordinate  system  remains  fixed.  The  parametric  formulation  of  the  curve  propagation  is 
in  a  Eagrangian  framework  in  which  the  coordinate  system  depends  on  the  parameterization  corre¬ 
sponding  to  the  curve’s  rest  position.  In  Eigure  2.1,  it  can  be  observed  that  the  level  set  formulation, 
where  a  level  curve  is  embedded  into  a  higher  dimensional  level  set  function  $,  may  be  more  advan- 
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tageous  in  handling  topological  changes  when  compared  to  parametric  evolution  methods  of  curves 
(also  called  marker  particle  methods  [126]).  This  is  due  to  the  fact  that  the  zero  level  set  of  $  need 
not  be  simply  connected,  may  split,  and  merge  during  the  evolution. 


Figure  2.1:  Depiction  of  a  level  set  formulation  versus  a  parametric  formulation  of  a  curve.  Left  column: 
(Eulerian  framework)  Evolution  of  a  level  set  function  ^{x,y,t)  as  a  graph  on  a  fixed  grid  is  shown  at  time 
instants  ti  and  t2 .  Note  that  topology  changes  are  handled  naturally;  Right  column:  (Lagrangian  framework) 
Evolution  of  zero  level  set  values  of  the  surface  via  marker  particle  methods,  i.e.  evolution  of  position  vector 
of  a  curve,  {X{p),y{p))  =  C  (p)  is  also  shown  at  time  instants  ti  and  f 2 ■  Notice  that  it  is  hard  to  manipulate 
topology  changes,  and  serious  bookkeeping  is  required. 

The  level  set  function  $  is  usually  picked  as  the  signed  distance  function  from  the  zero  level 
curve  which  is  the  contour  that  is  to  be  evolved.  A  fast  level  set  method  so  called  as  narrowband 
technique  [2]  developed  for  propagating  interfaces  is  used  in  the  implementation.  To  speed  up 
the  curve  evolution  algorithms  we  develop  in  this  thesis,  a  periodical  re-initialization  of  the  signed 
distance  function,  i.e.  the  level  set  function,  is  also  carried  out  by  the  technique  developed  by  [152]. 

In  the  level  set  equation,  (2.19),  the  speed  function  F  may  be  given  hy  F  =  a  +  j3K,  with  an 
advection  (first  term  on  the  right),  and  a  diffusion  term  or  a  curvature  term  (second  term  on  the 
right).  With  a  simplest  flow  when  F  =  1, 


(2.22) 


a  curve  propagates  in  its  inward  normal  direction  at  each  point  with  a  unit  speed,  and  it  develops 
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singularities,  i.e.  points  that  are  not  differentiable,  in  finite  time  [126].  Once  such  points  develop, 
the  normal  is  not  defined  af  fhose  poinfs,  and  fhe  propagafion  becomes  ambiguous.  Thus,  in  order 
fo  continue  fhe  evolufion,  a  weak  solution  is  required.  Nofe  fhaf  a  solution  is  a  weak  solution  of  a 
differential  equation  if  if  salisfies  an  infegral  formulalion  of  fhe  equafion.  This  in  furn  implies  fhaf 
fhe  weak  solution,  as  a  pofenfial  solufion,  may  nof  require  fhe  same  degree  of  differenliabilily,  and 
allows  for  more  general  solufions  [126].  A  way  fo  obfain  such  a  solution  is  provided  by  Sefhian, 
fhrough  fhe  notion  of  an  enfropy  condition,  also  called  Huygen’s  consfrucfion.  This  is  mofivafed  by 
an  analogy  of  a  curve  wifh  a  propagating  flame,  and  once  a  poinf  is  ignifed  by  fhe  expanding  flame 
(curve),  if  slays  burnl.  Thai  is,  information  once  losl,  can  nof  be  re-crealed  during  fhe  evolufion.  The 
parallelism  of  Ihis  developmenl  and  fhe  Iheory  of  viscosity  solutions  of  Hamillon-Jacobi  equalions, 
as  well  as  shocks  and  rarefaction  fans  formalions  in  hyperbolic  conservalion  laws  can  be  found 
in  [126].  To  solve  fhe  Eq.(2.22),  slable  numerical  schemes  are  developed,  which  selecl  fhe  correcl 
weak  solution  corresponding  fo  fhe  viscous  limils  of  fhe  associaled  curvalure-driven  equalions.  See 
[86,  88]  for  delailed  discussions  on  hyperbolic  conservation  laws.  Upwind  differencing  schemes  in 
computing  firsl  order  spatial  derivafives  which  use  values  upwind  of  fhe  direcfion  of  information 
propagation,  are  widely  employed  [126].  The  curve  evolution  equalions  we  develop  in  Ibis  Ihesis 
are  curvalure-driven  equations,  for  which  cenlral-differencing  schemes  for  fhe  spatial  derivatives 
and  forward-differencing  schemes  for  fhe  time  derivatives  are  adequate.  For  fhe  implemenlalion  of 
fhe  level  sef  melhod,  and  ils  re-inilializalion  slep,  however,  upwind  differencing  schemes  are  used. 

2.3  Overview  of  Segmentation  Methods 

The  problem  of  image  segmenlalion  refers  fo  fhe  partitioning  of  a  domain  of  a  given  image  I  info 
regions  such  fhaf  each  region  has  properties  disfincf  in  some  sense.  If  is  expecfed  fhaf  fhe  resulfing 
partitions  correspond  fo  meaningful  parls  of  objecls  in  an  image.  Image  segmenlalion  is  an  essential 
firsl  step  in  early  vision  and  provides  a  mechanism  for  an  automatic  analysis  of  image  conlenls. 
Some  imporlanl  applications  include  automatic  largel  recognition,  remote  sensing,  automatic  visual 
inspection  in  manufacluring  processes,  biomedical  image  analysis,  Iracking  objecls  in  motion,  and 
so  on.  In  Ihe  conlexl  of  remote  sensing  of  Ihe  earlh  for  inslance,  Ihe  image  would  be  partitioned  into 
regions  of  differenl  terrain  or  land  type. 

An  initial  lask  for  segmenlalion  is  to  determine  Ihe  fealures  which  delineate  and  reliably  dis¬ 
tinguish  differenl  regions.  Some  common  fealures  which  discriminate  Ihe  region  characteristics  are 
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intensity,  color,  and  texture.  Following  the  definition  of  features  for  the  segmentation  problem,  it 
is  necessary  to  select  a  “good”  criterion  for  capturing  and  evaluating  the  features  which  yield  a 
partition  of  the  image  domain  into  “different”  regions. 

Existing  image  segmentation  methods  may  broadly  be  classified  info  four  groups  [153]  as  de¬ 
scribed  in  fhe  nexf  subsecfions. 

2.3.1  Thresholding  and  Local  filtering  approaches 

Perhaps  fhe  earliesf  approaches  fo  image  segmenfafion  are  based  on  fhresholding  fechniques.  They 
rely  on  a  very  simple  concepf  which  is  fo  compare  each  pixel  value  {I{x,y),  {x,y)  G  fl)  wifh  a 
paramefer  (fhreshold)  and  decide  whefher  fhe  pixel  is  wifhin  fhe  region  or  nof.  The  value  of  fhe 
fhreshold  may  be  globally  or  locally  sef  [112, 134].  In  mosf  mefhods,  fhe  fhreshold  is  chosen  from 
fhe  infensify  histogram  of  eifher  fhe  whole  image  or  local  regions  of  fhe  image. 

Local  filfering  approaches  for  image  segmenfafion  are  based  on  defection  of  edges  which  cor¬ 
respond  fo  objecf  boundaries  or  fhe  boundaries  befween  image  regions.  The  celebrafed  early  work 
by  Marr  [96]  and  based  on  a  “primal  skefch”  concepf  enfailed  localizing  edges  in  fhe  image  for 
subsequenf  use  by  high  level  image  processing  steps.  Marr  and  Hildredfh  [97]  developed  an  edge 
defecfion  filfer  based  on  local  maxima  of  fhe  gradienf  magnifude.  When  fhe  firsf  derivative  achieves 
a  maximum,  fhe  second  derivafive  is  zero.  For  fhis  reason,  an  edge  defecfion  sfrafegy  is  fo  isolate 
zeros  of  fhe  second  derivatives  of  I.  The  differential  operafor  used  in  fhese  so-called  zero-crossing 
edge  defectors,  is  fhe  Laplacian.  In  order  fo  mitigate  fhe  increase  in  pixel  noise  due  fo  differenfi- 
afion,  fhe  image  is  pre-filfered  wifh  a  lowpass  filfer  such  as  a  Gaussian  kernel.  A  varianf  of  fhis 
idea  such  as  fhe  Canny  edge  defector  [22]  uses  fhe  zero-crossings  of  fhe  second  order  operafor  in 
Eq.(2.6)  insfead  of  fhe  Laplacian.  For  efficienf  implemenfafions.  Deriche  in  [40],  derives  exacf 
recursive  filfers,  faking  a  similar  analyfical  approach  fo  Canny. 

Segmenfafion  fechniques  in  fhis  group  make  use  of  local  informafion,  and  can  offen  be  imple¬ 
mented  as  a  convolution  of  fhe  given  image  wifh  fhe  impulse  response  of  a  local  filfer  for  efficiency 
sake.  They,  however,  rely  on  high  gradienf  values  in  fhe  image  fo  defecf  prominenf  boundaries 
befween  regions,  hence,  making  fhem  sensitive  fo  noise  and  affecting  fhe  confinuify  of  fhe  edge 
contours  as  a  resulf. 
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2.3.2  Region  growing  techniques 

Region  growing  techniques  partition  an  image  domain  into  k  disjoint  regions  Oi,  •  •  •  ,  Oh,  (i.e., 
O,  =  |jf=i  ^i)’  that  the  image  1°  is  homogeneous  in  some  sense  within  each  region.  Each 
region  Oi  is  a  connected  region.  Morel  and  Solimini  [100]  describe  a  general  multiscale  approach 
to  region  growing  as  follows:  (i)  Initialize  the  algorithm  with  the  finest  possible  segmentation  at 
a  small  scale  t,  i.e.  consider  each  pixel  as  a  separate  region;  (ii)  Merge  all  pairs  of  regions  whose 
“merging”  improves  the  segmentation;  (iii)  Iterate  (go  to  (ii))  by  increasing  the  scale  parameter. 
Choosing  the  criterion  to  perform  step  (ii)  results  in  different  algorithms. 

For  instance,  a  piecewise  (p-w)  constant  model  for  the  image  is  used  by  [78]  as  the  merging 
criterion.  To  a  region  Oi,  the  average  value  of  1°  over  that  region  is  assigned.  Given  average  values 
Ii  in  region  Oi  and  Ij  in  a  neighbor  region  Oj,  the  regions  Oi  and  Oj  are  merged  by  removing 
the  boundary  between  them  and  replacing  both  Ii  and  Ij  with  their  weighted  average  ■ 

Here  \Oi\  is  the  area  of  region  Oi .  The  algorithm  looks  for  a  decrease  of  global  energy  by  merging 
these  regions  and  by  updating  the  image  to  I.  The  global  energy  is  a  least  squares  criterion,  i.e. 
J  1 1 7  —  7°  I P  plus  a  penalty  on  the  total  length  of  the  boundaries  7.  However,  the  quadratic  penalty 
on  the  difference  between  the  estimate  7  and  the  initial  data  1°  is  not  well  suited  to  non-Gaussian 
noise,  such  as  speckle  noise  in  SAR  images. 

An  earlier  method  by  Pavlidis  in  [100]  also  starts  by  taking  all  pixels  as  regions,  and  merging 
every  pair  of  regions  Oi  and  Oj  such  that  variance  of  1°  over  Oi  |J  Oj  is  less  than  t  (scale). 

Variants  of  region  growing  or  merging  methods  may  yield  alternative  approaches,  e.g.  region 
splitting,  and  region  split-and-merge  methods  ( [100])  which  combines  the  two. 

Region  growing  methods  have  an  advantage  of  using  statistics  inside  regions,  they,  however, 
often  generate  irregular  boundaries  [153]. 

2.3.3  Active  contour  methods 

Active  contour  models  have  been  widely  used  in  image  segmentation  applications.  The  general  idea 
in  an  active  contour  framework  is  to  define  energy  funcfionals  whose  local  minima  comprise  a  sef 
of  solutions,  e.g.  fhe  boundaries  of  regions,  available  fo  higher  level  processes. 
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Snakes 


One  of  the  pioneering  works  in  this  field  is  due  to  Kass,  Witkin,  and  Terzopoulos  [70]  who  addressed 
the  problem  of  finding  salient  image  eontours  like  boundaries  of  objeets,  edges,  lines,  by  so-ealled 
“snakes”  algorithms.  They  aimed  at  having  the  snake  loek  onto  image  features  by  minimizing  an 
integral  measure  whieh  represents  the  snake’s  total  energy.  By  adding  suitable  energy  terms  to  the 
minimization,  it  is  possible  for  a  user  to  push  the  model  out  of  a  loeal  minimum  towards  the  desired 
solution.  Initially,  the  user  plaees  some  eontour  (snake)  near  an  image  strueture.  The  eonstraint 
forees  that  aet  on  a  snake  then  push  the  snake  towards  features  of  interest. 

Representing  the  position  of  the  eontour  parametrieally,  C  (p)  =  {X{p),y{p)),  where  p  G 
[0, 1],  Kass  et.al.  defined  the  snake’s  total  energy  funetional,  as 

1  .1 

.^image  {Cip))dp+  /  ^int(C(p))  T  .^con(C*  {p))dp.  (2.23) 

7o 

Here,  E'int  represents  the  internal  energy  of  the  snake  due  to  bending: 

=  (n;i(p)||Cp(p)||2  +  W2{p)\\Cpp{p)\?)l‘^. 

where  wi  and  W2  eontrol  the  “tension”  and  “rigidity”  of  the  snake  respeetively.  (Note  that  the 
subseripts  denote  derivatives  with  respeet  to  p,  and  1 1  •  1 1  denotes  the  standard  Euelidean  norm.) 
Their  basie  snake  model  is  a  spline  under  the  inlluenee  of  image  forees,  internal  eonstraints,  and 
other  general  eonstraint  forees,  E'con-  In  an  image  foree,  /g^  £'image(C'  {p))dp,  £'image(a:,  y)  is  a 
scalar  potential  field  defined  on  the  image  plane.  The  local  minima  of  J  E'image  attract  the  snake. 
For  example,  E'image  can  be  chosen  as  an  edge  functional  E'edge  =  “I  |V/(rE,  y)|p,  which  drives  the 
snake  to  contours  with  large  image  gradients,  i.e.  the  intensity  edges. 


E{C{p))=  [ 
Jo 


First  Connection  between  Curve  Evolutions  and  Active  Contours 


The  mathematical  foundation  of  another  class  of  geometric  active  contours  was  based  on  Euclidean 
curve  shortening.  Defining  the  length  functional  (see  arclength  definition  (2.9)) 
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then  differentiating  (taking  the  “first  variation”),  and  using  integration  by  parts,  we  have 


where  k  is  the  eurvature,  N  is  the  inward  unit  normal,  and  =  ds.  Thus  the  direetion  in 

whieh  L{t)  is  deereasing  most  rapidly  is  aehieved  when 


dC 


kN 


whieh  defines  a  gradient  flow.  (For  fhe  derivafion  of  fhis  eurve  lengfh  shorfening  flow,  see  [148].) 
A  new  aefive  eonfour  paradigm  was  proposed  in  [25, 72, 148],  by  ehanging  fhis  ordinary  Euelidean 
are-lengfh  funelion  along  a  eurve  C  (p)  wifh  paramefer  p  given  by 


ds  =  \\Cp\\dp  =  {x^  +  ylfi‘^dp 


fo 

dS(f,  =  (f)ds  =  {Xp  +  yp)^^‘^(f)dp 


where  y)  is  a  positive  differentiable  funelion.  Thus  a  new  melrie  is  defined  wifh  whieh  a  new 
gradienl  flow  is  lo  be  found.  The  gradienf  flow  for  fhe  eurve  shortening  relative  fo  fhe  new  melrie 
dS(j,  will  be  eompuled,  Ihen  fhe  firsl  variation  of 


will  lead  lo 


dC 


(</)K  -  V</)  •  iV  )iV . 


The  melrie  ds^  has  fhe  properly  lhal  il  beeomes  small  where  cj)  is  small  and  viee  versa.  Thus  al  sueh 
poinls,  lenglhs  deerease,  and  less  energy  is  needed  in  order  lo  move.  If  one  wanls  Ihe  eonlour  loek 
onto  edges  of  an  image,  Ihen  il  is  reasonable  to  eonslruel  a  weighl  whieh  is  almosl  zero  near  edges, 
and  almost  1  when  it  is  far  from  the  edges.  Sinee  ||V/||  is  a  loeal  indieator  of  strength  of  edges  in 
an  image,  (j)  is  ehosen  as 

^  =  TTW^\- 

These  geometrieal  snakes  are  very  loeal  models,  and  are  only  sensitive  to  data  near  the  eurve. 
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Figure  2.2:  A  piece-wise  constant  image  consisting  of  two  regions. 

The  initial  contour  has  hence  to  be  reasonably  close  to  boundaries.  The  local  property  also  causes 
these  models  to  be  very  sensitive  to  noise. 

Region-based  Active  Contours 

Region-based  active  contours  were  proposed  to  overcome  the  problems  with  the  geometric  or 
geodesic  [25, 72, 148]  active  contours  by  using  both  local  and  global  information.  The  main  idea  is 
to  assume  that  the  image  consists  of  a  finite  number  of  regions,  that  are  characterized  by  a  prede¬ 
termined  set  of  features  or  statistics  such  as  means,  variances,  textures  etc.  These  features/statistics 
can  be  estimated  from  the  image  data,  hence  an  energy  functional  may  be  constructed  to  pull  these 
statistics  apart,  i.e.  maximize  the  distance  between  them  in  order  to  separate  the  corresponding 
regions.  One  advantage  over  the  models  in  the  previous  subsection  is  that  there  is  no  need  to  cal¬ 
culate  gradients  of  the  image  which  are  usually  very  sensitive  to  noise.  Region  based  flows  are 
therefore  much  more  robust  to  noise,  at  a  cost  of  additional  imposed  assumptions  on  the  images, 
and  additional  computations. 

One  such  assumption  is  one’s  ability  to  approximate  an  image  by  constants,  i.e.  assuming  that 
the  image  consists  of  piecewise  constant  regions.  Let  an  image  consist  of  only  two  regions,  a 
foreground  Rj,  and  a  background  0,\Rf,  and  let  these  be  approximated  by  constants.  Fix  C ,  an 
arbitrary  closed  curve  in  the  image  domain  Q,  and  therefore  0,  is  partitioned  into  R  and  i?®,  regions 
inside  and  outside  C  respectively  (Figure  2.2).  An  energy  functional  to  separate  the  means  of  the 
two  regions  R  and  R^,  say  u  and  v  respectively,  is  given  by  Yezzi,  Tsai,  and  Willsky  [149] 

Eytw{C  )  =  --(t(  -  (2.24) 

A  gradient-descent  flow  on  the  above  energy  functional  will  maximize  the  Euclidean  distance  be¬ 
tween  the  mean  of  region  R  and  the  mean  of  region  Q\R.  The  energy  in  (2.24)  hence  will  be 
minimum  when  the  curve  C  locks  onto  the  boundary  between  the  foreground  (target  object),  and 
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the  background. 

A  similar  energy  functional  is  given  by  Chan  and  Vese  [29] 

Ecv(C)=  [  [  (I-ufdxdy+  [  [  {I-vfdxdy.  (2.25) 

J  Jr  J  Jn\R 

The  resulting  flow  in  the  direction  of  the  gradient  descent  will  automatically  move  the  contour 
towards  the  boundary  dRj  of  the  target  shape,  to  minimize  this  fitting  energy.  For  instance,  an 
initial  contour  encompassing  the  region  Rj  will  flow  inward  towards  the  boundary;  while  a  contour 
inside  Rf  will  flow  outward;  and  a  contour  which  overlaps  Rj  will  flow  in  both  directions  towards 
the  boundary.  This  makes  the  initial  placement  of  the  contour  less  restrictive. 

The  data  term,  resulting  from  this  formulation,  was,  however,  shown  to  possibly  require  addi¬ 
tional  regularization  depending  on  the  particular  types  of  noise,  e.g.  for  salt  and  pepper  type  noise, 
the  contour  may  weave  around  noisy  regions  and  result  in  erroneous  regions.  A  regularizing  term 
(penalty  on  the  length  of  the  curve)  is  added  to  yield  for  instance  the  following  energy  functional 


E(C )  =  F/ytw  +  cx. 


(2.26) 


where  ds  is  the  total  arclength  of  the  curve,  and  tv  is  a  parameter  which  determines  the  amount 
of  the  desired  regularization.  The  aim  is  to  thus  prevent  the  length  of  the  curve  from  getting  imprac- 
tically  long  and  producing  smooth  boundaries.  The  corresponding  gradient  descent  flow  for  E'ytw, 
can  be  shown  to  be 


dC  ,  ,,  I  -  u  I  -  V 

—  =  (u-v){  +j - ——)N-aKN, 

c)t  J  ji  QiX(ty 

=  fiN  -  anN  .  (2.27) 

(A  very  similar  derivation  for  the  first  term  (data  term)  is  deferred  to  Chapter  4).  The  regularizing 
term  (second  term  above)  may  also  be  viewed  as  a  shape  prior  which  is  especially  very  strong  at 
contour  points  with  high  curvature.  Therefore,  there  is  a  trade-off  between  the  data-driven  term  and 
the  regularizing  term.  The  data-driven  term,  however,  usually  dominates  the  flow. 

2.3.4  Global  optimization  approaches  based  on  energy  functionals 

The  goal  of  most  active  contour  algorithms  is  to  extract  the  boundaries  of  homogeneous  regions 
within  an  image,  whereas  the  goal  of  most  nonlinear  diffusion  algorithms  is  to  smooth  an  image 
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within  homogeneous  regions  but  not  aeross  the  boundaries  of  sueh  regions.  A  well-known  math- 
ematieal  model  proposed  by  Mumford  and  Shah  [102]  simultaneously  addresses  both  goals.  They 
develop  an  energy  funetional  whieh  approximates  an  image  by  smooth  funetions  in  eaeh  region 
instead  of  eonstant  ones.  The  Mumford-Shah  (M-S)  funetional  is  given  by 

E{CjR,f„,)  =  JI(,fx-lfdxdy  +  JI{fE.-lfdxdy 

where  C  is  the  elosed,  smooth  segmenting  eurve,  I  is  the  observed  image  data,  is  the  smooth 
funetion  inside  the  eurve,  f^c  is  the  smooth  funetion  outside  the  eurve.  Minimizing  Eq.  (2.28)  then 
eorresponds  to  finding  estimates  /r  and  f^c  in  regions  R  and  respeetively. 

The  first  two  terms  in  the  M-S  funetional  are  the  data-fidelity  terms  (like  the  measurement/ 
observation  model),  the  seeond  two  terms  are  the  smoothness  terms  in  the  given  regions  (like  a 
prior  model  for  /  given  C ).  The  last  term  is  a  prior  model  for  C  whieh  penalizes  its  are  length. 
The  M-S  funetional,  henee,  eaptures  the  desired  properties  of  segmentation  and  reeonstruetion  by 
pieeewise  smooth  funetions  as  opposed  to  p-w  eonstant  models  in  Eq.(2.24),  and  Eq.(2.25). 

There  are  numerous  other  algorithms  based  on  minimizing  different  eriteria  sueh  as  Minimum 
Deseription  Eength  (MDE)  eriteria  [87],  Bayesian  eriteria  [17].  The  problem  of  minimizing  energy 
funetionals  sueh  as  Eq.(2.28)  or  other  energy  funetionals  (MDE  or  Bayes  eriteria)  is  prineipally 
eomputational  (e.g.  simulated  annealing  [51]).  To  overeome  this  diffieulty,  algorithms  sueh  as 
mean  field  annealing  [14,  61, 132],  graduafed  non-eonvexify  [17]  have  also  been  explored. 
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Chapter  3 


Stochastic  Differential  Equations  and 
Geometric  Flows 


In  this  chapter,  we  present  a  new  class  of  curve  evolution  equations  for  smoothing  of  structures 
along  the  known  orientation  of  salient  lines  in  curves  and  images  while  preserving  their  important 
features.  The  contents  of  this  chapter  are  outlined  as  follows.  After  an  introduction  in  Section  3.1, 
we  briefly  review  and  recap,  as  was  descrihed  in  Chapter  2.2,  some  theoretical  concepts  associated 
with  the  curve  shortening  flow,  including  its  connection  to  a  nonlinear,  directional  diffusion  equa¬ 
tion  in  which  image  values  diffuse  locally  only  along  the  directions  of  its  edges  in  Section  3.2.  In 
Section  3.3,  we  provide  a  stochastic  equivalent  equation  which  in  turn  unveils  a  new  shape/feature- 
driven  flow  descrihed  in  detail  in  Section  3.4,  which  we  also  believe  could  offer  a  variety  of  appli¬ 
cations  outside  the  recognition  and  classification  problems.  We  conclude  with  some  illustrating  and 
substantiating  examples  in  Section  3.5,  and  conclusions  in  Section  3.6. 

3.1  Introduction 

In  recent  years  curve  evolution  has  emerged  as  an  important  application  of  partial  differential  equa¬ 
tions  (PDE’s)  in  image  processing,  computer  vision,  and  computer  graphics.  Curve  evolution  tech¬ 
niques  have  been  applied  not  only  to  individual  curves,  for  applications  such  as  edge-detection, 
skeletonization,  and  shape  analysis,  but  have  also  been  considered  for  their  simultaneous  action  on 
the  level  sets  of  an  image  in  a  number  of  geometrically  based  anisotropic  smoothing  algorithms. 
Osher  and  Sethian  [111,  126]  extended  this  latter  perspective  to  the  treatment  of  individual  curves 
through  a  set  of  algorithms,  known  as  level  set  methods,  which  enable  the  implementation  of  curve 
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and  surface  evolution  on  a  fixed  grid.  These  techniques  have  aided  a  number  of  researchers  in  push¬ 
ing  the  application  of  curve  evolution  to  new  limits  by  providing  a  simple  framework  for  treating 
certain  types  of  singularities  such  as  shocks  and  topological  transitions  [109,  111]. 

Much  of  the  research  in  curve  evolution  theory  has  centered  around  the  so  called  geometric  heat 
equation  [54]  in  which  a  curve  is  evolved  along  the  normal  direction  in  proportion  to  its  signed 
curvature.  This  flow  is  well  known  for  its  smoothing  properties  [74-76]  and  the  fact  that  it  cor¬ 
responds  to  the  gradient  evolution  for  arclength  (thereby  earning  the  name  curve  shortening  flow). 
Because  curvature  is  a  purely  geometric  quantity  (invariant  to  rotation  and  translation),  curvature- 
based  motion  gives  rise  to  a  Euclidean  invariant  scale  space  [5, 6, 147],  allowing  one  to  trace  fea¬ 
tures  in  a  curve  from  finer  to  coarser  scales  as  the  evolution  proceeds.  An  affine  invariant  scale  space 
can  be  obtained  from  a  related  curvature  flow  which  depends  upon  the  cube  root  of  the  curvature 
(see  [5, 124, 125]). 

When  applied  to  the  level  sets  of  an  image,  these  flows  have  a  powerful  denoising  effect  when 
run  for  a  short  amount  of  time.  If  run  for  too  long,  however,  even  large  scale  features  will  be  de¬ 
stroyed.  The  reason  stems  from  the  fact  that  as  the  geometric  heat  flow  shrinks  any  closed  curve,  the 
curve  becomes  more  and  more  circular  (elliptical  in  the  case  of  the  affine  flow)  and  will  eventually 
collapse  into  a  single  point  [54].  It  is  therefore  not  always  possible  to  preserve  desired  features  in 
the  shapes  of  objects  (corners  for  example)  if  too  much  evolution  is  required  to  remove  a  significant 
level  of  noise.  Furthermore,  it  is  not  well  understood  how  these  curvature-based  filters  are  affected 
by  different  noise  distributions  and  when  this  sort  of  problem  may  occur. 

To  the  best  of  our  knowledge,  and  aside  from  [80,  81],  nonlinear  diffusion  in  the  previous  litera¬ 
ture  was  discussed  from  a  purely  deterministic  perspective.  In  this  chapter,  we  provide  a  stochastic 
formulation  of  the  geometric  heat  equation  and  use  the  resulting  insights  to  develop  a  new  class  of 
curvature-based  flows  and  anisotropic  diffusion  filters  which  preserve  desired  features  in  the  shape 
of  an  object.  Under  these  new  flows,  evolving  curves  take  the  limiting  form  of  a  polygon  (see  [20] 
for  evolutions  of  polygons  related  to  the  geometric  and  affine  geometric  heat  flows,  and  [144]  for 
evolutions  of  polygons  globally  through  an  electric  field  concept).  The  resulting  diffusion  models 
may  therefore  be  applied  for  much  longer  periods  of  time  without  distorting  the  shapes  of  polygonal 
objects  in  the  image,  thereby  mitigating  the  tradeoff  between  noise  removal  and  shape  distortion. 

Polygonal  structures  are  ubiquitous  in  images  of  man-made  objects  (buildings,  roads,  vehicles, 
and  so  on),  which  contain  many  straight  lines,  often  oriented  in  particular  directions  (e.g.  horizontal 
and  vertical),  that  come  together  to  form  sharp  corners.  The  ability  to  preserve  such  distinctive 
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features  is  not  only  desirable  when  filtering  an  image  which  contains  these  types  of  shapes,  but 
is  also  important  when  applying  low  level  smoothing  to  an  extracted  shape  since  such  features 
constitute  important  and  powerful  cues  for  recognizing  objects  in  higher  level  vision  algorithms. 
We  will  present  both  applications  in  this  chapter.  From  a  dual  perspective  to  our  contour-based 
approach  to  shape  representation,  skeletonization  approaches  may  also  allow  shape  analysis  without 
displacement  of  corners  [18, 35, 75, 107, 118, 129]. 

In  this  chapter,  we  develop  a  new  class  of  curve  evolutions,  which  are  obtained  by  a  modification 
of  the  geometric  heat  equation.  Given  an  initial  shape  in  the  form  of  a  continuous  curve,  the  class  of 
curve  evolution  equations  we  will  obtain,  deform  it  into  a  pre-specified  final  polygonal  shape.  The 
problem  of  deforming  an  inpuf  shape  info  a  differenf  form  has  been  of  inferesf  in  various  fields  such 
as  compufer  graphics  [52]. 

3.2  Background 

The  geomefric  heal  equation,  infroduced  in  Chapter  2,  may  be  oblained  as  can  be  recalled  in  a 
geomefric  manner  by  specifying  fhe  evolufion  of  each  level  curve  in  fhe  image.  Lef  C  denote  a 
parficular  iso-infensify  confour  which  we  will  deform  over  time  via  fhe  following  flow, 

Ct  =  Css  =  nN  (3.1) 

where  s  denofes  fhe  arclengfh  paramefer,  k  fhe  Euclidean  curvalure,  and  N  fhe  inward  unif  normal. 
Equation  (3.1),  referred  fo  as  fhe  Geometric  Heat  Equation  (GHE),  is  well  known  for  ifs  smoofhing 
properties.  If  has  been  shown  by  Grayson  [54]  fhaf  any  closed,  embedded  curve  evolving  according 
fo  (3. 1)  will  convexify  and  smoofhly  shrink  lo  a  single  poinf  in  finite  time,  becoming  more  and  more 
circular  along  fhe  way.  If  we  apply  fhe  geomefric  heal  flow  lo  every  single  level  curve  in  fhe  image 
we  oblain  fhe  same  anisolropic  diffusion  equation  lhal  we  derived  in  Chapter  2.  To  see  Ihis,  note 
lhal  al  lime  t  each  level  curve  C  ^  (where  fhe  index  k  distinguishes  one  level  curve  from  anolher) 
is  implicifly  described  by  u{t,  x,  y)  =  where  denofes  a  parficular  infensily  in  fhe  image.  Eel 
us  choose  a  parameterization  of  so  lhal  y{t,p))  for  p  €  [0, 1]  and  for 

all  f  >  0.  We  may  Ihen  write  u{t,C^{t,p))  =  u(t,  X{t,p),y{t,p))  =  u^.  Differentiating  Ihis 
expression  wilh  respecf  lo  t  yields 

ut  +  Vu  ■  C  t  =  ut  +  Vu  ■  {kN  )  =  0. 
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Note  that  the  inward  unit  normal  and  the  eurvature  of  eaeh  level  eurve  ean  be  expressed  as  AT  = 
—  and  K  =  V  •  ^  ^ .  This  allows  us  to  rewrite  the  above  equation  eompletely  in  terms  of  u 

and  its  derivatives, 


ut  =  V 


f  Vu 

■  IvIlVull 


llVull  = 


u: 


I'^xx  ‘^'^x'^y'^xy  T 


x'^yy 


Ux  T  Uy 


(3.2) 


giving  us  a  PDE  whieh  is  identieal  to  (2.8). 

Equation  (3.2)  is  also  referred  to  as  the  geometrie  heat  equation  sinee  it  eomes  from  applying 
the  previous  geometrie  heat  equation  (3.1)  to  eaeh  level  eurve  of  an  image  u.  This  double  meaning 
of  the  term  geometric  heat  equation  is  disambiguated  by  the  eontext  in  whieh  the  flow  is  applied 
(i.e.  either  to  an  image  or  to  a  eurve).  In  this  ehapter,  we  will  be  interested  in  both  eases  and 
will  present  direetional  generalizations  of  the  geometrie  heat  flow  whieh  are  designed  to  preserve 
eertain  types  of  features  either  in  a  eurve  or  in  an  image.  We  first  however,  reformulate  the  geometrie 
heat  flow  from  a  stoehastie  point  of  view,  giving  new  insights  into  the  nature  and  behavior  of  this 
nonlinear  diffusion  model.  It  was  preeisely  these  insights  that  led  us  to  the  generalizations  presented 
in  Seetion  3.4. 


Remark:  We  note  that,  in  general,  the  popular  anisotropie  diffusion  models  of  Perona  and  Malik 
in  [117],  are  not  related  to  eurve  evolution  theory  and  are  only  intended  for  images,  not  eurves 
(unless  the  eurve  has  the  form  of  a  graph).  As  sueh,  we  will  not  attempt  to  relate  the  eurve  evolution 
models  developed  in  this  thesis  to  Perona-Malik  models  whieh  represent  a  different  perspeetive  on 
the  subjeet  of  nonlinear  diffusion. 


3.3  Stochastic  Formulation  of  a  Geometric  Heat  Equation 

3.3.1  Introduction  to  Ito  Diffusion 

The  diffusion  of  a  partiele  is  usually  well  modeled  by  a  Stoehastie  Differential  Equation  (SDE) 
whieh,  in  turn,  represents  the  underlying  mieroseopie  proeess  of  an  evolution  of  a  pixel  or  a  point. 
The  dynamies  of  this  evolution  at  a  maeroseopie  level  are  eaptured  by  a  PDE,  also  referred  to  as 
a  generator  (infinitesimal)  of  the  diffusion  [80,  81, 108].  Suppose  we  want  to  deseribe  the  motion 
of  a  small  partiele  suspended  in  a  moving  liquid,  subjeet  to  random  moleeular  bombardments.  If 
5  (f,  a: )  G  M”  is  the  veloeity  of  the  fluid  at  a  point  a:  G  M”  and  time  t  G  M+,  then  a  widely  used 
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mathematical  model  for  the  position  X  (t)  of  the  particle  at  time  t  is  an  SDE  of  the  form 


dX  (t)  =  b{t,X  {t))dt  +  (T{t,X  {t))dB  (t), 


(3.3) 


where  X  (t)  is  an  n, -dimensional  stochastic  process,  a  {t,x)  G  ,  and  B  (t)  is  an  tt?, -dimensional 
Brownian  motion.  6  (•,  •)  is  called  the  drift  coefficient,  and  cr  (•,  •)  is  called  the  diffusion  coefficient. 
The  first  term  in  this  equation  corresponds  to  a  non-random/deterministic  motion,  whereas  the  sec¬ 
ond  term  models  randomness  or  noise  in  the  motion. 

The  solution  of  such  an  SDE  may  be  thought  of  as  a  mathematical  description  of  the  motion  of  a 
small  particle  in  a  moving  fluid,  and  such  stochastic  processes  are  called  (Ito)  diffusions  [108].  Eor 
many  applications,  a  second  order  partial  differential  operator  A  can  be  associated  to  an  Ito  diffusion 
X  (t)  given  by  Eq.  (3.3).  The  basic  connection  between  A  and  X  (t)  is  that  A  is  the  generator  of 
the  process  X  (t).  \fw{x  )  G  Cq  (M”),  (i.e.,  it  is  continuous  with  continuous  derivatives  up  to  order 
2,  and  has  a  compact  support),  then  A  is  given  in  the  form 


A  w 


1 

2 


d‘^w  ^  dw 

i 


(3.4) 


In  conjunction  with  this,  the  so-called  Kolmogorov’s  backward  equation  [108],  gives  a  probabilistic 
solution  to  linear  partial  differential  equations.  Kolmogorov’s  theorem  states  that  given  X  (t)  = 
(2l  (f),  2l  (t)),  where  [•]  is  the  expectation  operator  with  respect  to  the  probability  law  of 

X  (t)  starting  at  the  point  x  ,  and  defining  7(f,  x  )  =  [f{X  (i))]>  then  there  exists  an  operator 
A  such  that 

d'y  o 

=  ,4-7,  f  >  0,  X  G 

ot 

7(0,  a:)  =  f{x),  ojGM^. 


SDEs  and  stochastic  processes,  most  commonly  the  Brownian  motion,  have  previously  been 
used  in  curve  and  image  analysis.  Mum  ford  [101]  used  it  to  model  completion  curves  of  occluded 
edges,  the  so-called  elastica.  By  taking  the  curvature  function  (of  arc  length)  as  a  Gaussian  process, 
and  the  tangent  direction  on  the  curve  then  as  a  Brownian  motion,  he  derived  the  probability  of  the 
curves  that  link  occluded  edges.  Eor  a  more  general  situation,  e.g.  curves  in  ,  Mumford  used  other 
sorts  of  stochastic  processes  such  as  an  Uhlenbeck  process  to  find  the  elastica.  Williams  and  Jacobs 
[146],  later  in  their  “Stochastic  Completion  Eields”  work,  define  the  same  SDE  as  Mumford’s,  for 
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a  particle’s  position  and  the  orientation,  and  through  this  model  of  diffusion  incorporate  the  prior 
assumption  that  the  maximum  likelihood  path  followed  by  a  particle  between  two  positions  and 
directions  is  a  curve  of  least  energy,  and  solve  it  by  a  discrete  formulation.  Similarly,  a  Kalman  filter 
which  produces  estimates  of  a  system  as  it  evolves  in  time  and  affected  by  noise,  (which  is  indeed  an 
SDE  written  for  the  system  and  its  observations),  was  used  in  [37]  for  grouping  of  contour  segments. 
Our  use  of  SDEs  is  along  a  different  line  of  thought  in  that  our  inspiration  starts  with  a  desired  effect 
of  a  nonlinear  filter.  Specifically,  fhe  fheory  of  SDEs  provides  us  wifh  a  microscopical  inferprefafion 
of  fhe  well-sfudied  geomefric  heal  equalion,  and  leads  lo  a  new  macroscopic  description  of  Ibis 
equation  which  in  furn  is  used  lo  develop  a  new  class  of  curve  evolutions  or  fillers. 

3.3.2  Stochastic  Formulation  of  the  Geometric  Heat  Equation 

Eel  us  denofe  by  0{t,x)  fhe  angle  belween  fhe  oulward  normal  lo  Ihe  curve  and  Ihe  x-axis  al 
each  spatial  poinl  x  =  {x,y).  The  oulward  unil  normal  N  can  Ihen  be  expressed  in  terms  of 
Ihe  angle  0  as  N  =  {cos{0{t,  x  ),  sin0(l,  x  )),  which  is  re-wrillen  in  terms  of  u(-)  as  N  {t,x)  = 
{Ux{t,  X  ),Uy{t,  X  ))  I Ux{t,xf  +  Uy{t,xf .  Ilfollows,  9{ux{t,  x),Uy{t,x))  =  tan~^(“^|*’^|). 
Using  Ihese  equations,  and  defining  an  operator  .Aghe  of  Ihe  form 

Aaw.u{t,x)  =  sir?  e{ux{t,x),uy{t,x))uxx{t,x) 

—  2  sin6{ux{t,  X  ),Uy{t,  x  ))  cos  6{ux{t,  x  ),Uy{t,  x  ))  Uxy{t,  x  ) 

+  COS^  6{Ux{t,  X  ),Uy{t,  X  ))  Uyy{t,  X  ),  (3.5) 

Ihe  geomelric  heal  equation  (3.2)  can  be  re-wrillen  as 

u{0,x)  =  Uo{x), 

Ut{t,x)  =  AaHEu{t,x),  (3.6) 

where  Uo{x  )  is  Ihe  initial  level  sel  function. 

In  lighl  of  Ihe  foregoing  developmenl,  a  nalural  question  which  arises  is:  given  a  PDE  which 
governs  a  curve  shortening  flow,  can  we  obtain  a  corresponding  SDE  associated  with  the  underlying 
diffusion  ? 

The  nonlinearily  of  GHE  presenls  a  signilicanl  challenge  to  find  a  global  Ilo  diffusion  which 
explains  Ihe  overall  microscopical  behavior  of  Ihe  system.  Our  approach  here  for  solving  such  a 
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nonlinear  problem  is,  to  explore  the  short-time  behavior  by  linearizing  around  a  known  (nominal) 
solution.  The  perturbation  equations  so  obtained  will  be  linear  and  henee  an  approximate  solution 
to  the  nonlinear  problem  ean  be  obtained  as  the  nominal  value  plus  the  perturbation  term.  Let  us 
denote  by  u^{t,  x  )  the  solution  to  Eq  (3.6): 


du^ 


sin^ 


U^y  +  COS^ 


Uyy, 


and  if  we  write  u{t,x)  as 

u{t,  X  )  =  vP{t,  x  )  +  e  u{t,  X  ), 

1  u^it  X ) 

and  define  the  eorresponding  nominal  angle  0^{t,x  )  =  tan~^(^j(^^’^^),  we  get  a  linearized  version 
of  the  geometrie  heat  equation  around  a  nominal  value: 


^  ^  «  AGHEMU{t,x  ) 

ot 

=  sva^{9‘^{x  ))  Uxx{t,  X  )  —  sin(20”(a: ))  Uxy(t,  x  )  +  cos^(6^(x  ))  Uyy(t,  x  ) 

+  c{x){-Uy{x)  Ux{t,x)  +Ux{x)  Uy{t,x)),  (3.7) 

where  c(a:)  =  ^uUx)rliu”ix)r  -  Uyy{x))  -  cos(2e^(x))2u^y(x)]  (see 

Appendix  A  for  details  of  this  derivation). 

In  light  of  this,  we  ean  proeeed  to  state  the  following: 


Proposition  1  The  right  hand  side  of  the  linear  PDE  in  Eq.  (3.7)  is  the  generator  of  the  following 
Ito  diffusion  satisfying  the  SDE 


I  dX(i)(f) 

dX  (2)  (t) 


c{X  (t)) 


-u-iX  (t)) 

Kix  it)) 


dt-\~  s/2 


-sine^iX  it)) 
cos  e^ix  it)) 


dBit). 


(3.8) 
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Proof  1  The  operator  AoHEUn  in  Eq.  (3.7)  is  first  re-written  as, 


AoHEUn  =  b^{X)-V  +  ^(t{X)ct^{X)  OH 

V  <(^)  /  \fy 

(  sin2  e^{X)  -  sin  (X  )  cos  (X  ) 

y  -sinr(X)cosr(X)  cos2r(X) 


©H, 


where  H  is  a  Hessian  operator  and  ©  is  a  Hadamard  product.  The  factorization  of  la  leads  to 


a{X)  =  s/2 


-sin0”(X) 
cos  e^{x) 


and  by  identification,  b  (X  ) 


/  -n;(X) 

V  <(^) 


Given  the  functions  b  (X  (t)),  and  a  (X  (t)),  we  come  up  with  a  pair  of  processes  (X  (i),  B{t)) 
such  that  the  SDE  in  Eq.  (3.8)  holds.  In  this  case,  the  solution  X  (t)  is  called  a  weak  solution,  as 
it  does  not  specify  beforehand  the  explicit  representation  of  the  white  noise,  i.e.  the  version  B(t)  of 
the  Brownian  motion  is  not  given  in  advance. 


Both  the  drift  and  diffusion  coefficient  vectors  of  this  SDE  are  in  the  tangent  direction  of  our 
level  curves,  which  helps  us  interpret  it  as  a  1 -dimensional  Ito  diffusion  on  the  instantaneous  tangent 
direction  T  {u^{t),Uy{t)).  A  differentiability  assumption  on  u{t,x) 

u{t  +  St,x)  —  u{t,x)  du{t,x) 
hm  ^ ^  ;  «  AoHEMUit,  X 

dt^o  5t  dt 

is  sufficient  for  a  short-time  existence  of  the  linearized  PDE  version  of  the  nonlinear  geometric  heat 
equation. 

Using  Kolmogorov’s  theorem  cited  in  Section  3.3.1,  and  assuming  that  u{t,  x  )  and  its  deriva¬ 
tives  are  “sufficiently  regular”  (Eipschitz  properties),  starting  at  each  time  t,  the  diffusion  X  (t)  in 
Eq.  (3.8)  is  constructed  for  each  time  interval  (t  —  5t,t),  and  may  be  used  to  write  a  Backward 
Kolmogorov  Equation, 


u{t  —  5t,x)  =  E{u{t,  X  {t))/X  (t  —  5t)  =  X  }, 
as  a  mean  value  around  each  pixel  dictated  by  the  motion  of  the  constructed  diffusion  process 
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X  (t).  This  equation  can  also  be  written  in  forward  time  (since  in  the  small  time  step  5t,  the  ap¬ 
proximate  constant-coefficient  PDE  gives  rise  to  a  time-homogeneous  diffusion  X  (t)  with  b{X  (t)) 
and  <T  (X  (t))  to  give  way  to  an  averaging  process  in  the  tangent  direction  of  a  level  curve  in  the 
course  of  a  forward  evolution  (i.e.,  estimate  the  new  pixel  value  at  time  f  as  a  mean  value  of  two 
neighboring  pixel  values  on  the  tangent  at  time  (t  —  5t)  (See  Fig.  3.1) ). 


Figure  3.1:  Illustration  of  averaging  behaviour.  Points  of  the  zero-level  set,  i.e.  initial  contour 
(T’(f),  y{t)),  at  time  t,  is  shown  on  the  left.  Those  points  whose  sample  realizations  result  in  an 
average  value  of  zero  at  time  t+5t  {u{t+5t,  x  )  =  (f))]  =  0)  form  the  new  contour(T’  {t+ 

St),  y{t  +  St))  (on  the  right). 


This  also  leads  us  to  infer  that  locally,  we  can  write  a  valid  diffusion  for  each  time  interval 
{t  —  St,  t), 

dX  (t)  =  T^(X  (t))  (c^(X  (t))  dt  +  V2  dB{t))  (3.9) 


where  T”(X  (t))  denotes  the  known  tangent  vector  at  time  t,  and  c”(X  (t))  is  the  known  drift 
coefficient  at  time  t,  which  is  c”  =  ^  ==  [sin(20”)(r/”  —  uf.,,)  —  cos(20”)2r/”  1  .  Apart 

from  a  drift  on  T ,  i.e.. 


dX  {t)  =  y/2T^{X  (t))  dB{t), 


(3.10) 


the  underlying  particle  motion  can  be  interpreted  as  a  Brownian  Motion  (BM)  on  a  local  frame  in 
the  direction  of  the  tangent.  Since  BM  is  an  averaging  process,  this  SDE  explains  that  the  geometric 
heat  equation  smooths  iso-intensity  contours  maximally.  On  a  discrete  lattice  Brownian  motion  is 
captured  by  a  random  walk  with  equally  likely  (i.e.,  prob.  1/2)  displacements  to  and  uTp.  The 
latters  are  obtained  by  a  bilinear  interpolation  around  u  (in  x  and  y  direction  and  along  the  tangent. 
See  Fig.  3.2). 

T  0.5 


Figure  3.2:  Symmetric  random  walk  on  the  tangent  direction,  and  corresponding  interpolation  on 
square  grid. 
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As  a  result,  we  can  write  such  an  equation  as 


du 

'm 


+ 


(3.11) 


where  At  is  the  Laplacian  operator  in  a  tangent  direction  T .  Summing  up,  microscopic  dynamics 
of  the  system  captured  by  the  local  diffusion  X  (t)  lead  to  a  new  macroscopic  description  of  the 
scenario,  i.e.,  the  random  walk  obtained  in  terms  of  the  macroscopic  variable  u{t,x). 

For  simulation  purposes,  we  use  a  level  set  methodology  [111],  which  in  an  Eulerian  framework 
has  an  advantage  of  naturally  handling  topological  changes  on  the  level  set  function.  A  simulation 
example  where  a  “T”  shape  is  evolved  via  a  generator  of  a  random  walk  in  a  tangent  direction, 
Eq  (3.11),  is  shown  in  Eig.  3.3. 


Eigure  3.3:  Generator  of  symmetric  random  walk  on  the  tangent  direction  implemented  on  the  level 
set  function  u{x,  y)  for  a  T  shape.  The  tangent  direction  is  estimated  directly  from  the  level  set  set 
function  :  Ot  =  tan“^  level  set  function  is  on  a  250  x  250  grid,  8t  =  0.25. 

Practical  equivalence  of  GHE  and  random  walk  on  the  tangent  direction  is  tested  by  several 
shapes.  Another  illustrating  simulation  is  shown  in  Eigure  3.4. 

Our  neglecting  the  drift  led  to  an  unbiased  random-walk  on  the  tangential  direction  and  is  val¬ 
idated  by  the  simulation  examples  presented  above,  as  the  generator  of  symmetric  random-walk 
implementation  results  are  in  agreement  with  the  geometric  heat  equation  implementation.  Theo¬ 
retically,  a  stronger  validation  is  due  to  Girsanov  theorem  (see  [108]),  which  says  that  if  we  change 
the  drift  coefficient  of  a  given  Ito  process,  then  the  law  of  this  process  does  not  change  dramatically, 
indeed,  the  trajectories  of  the  process  (distribution)  change  via  the  measure  change  on  the  trajecto¬ 
ries.  This  theorem  involving  a  change  of  measure  provides  us  with  a  means  of  changing  the  mean  of 
the  process  X  (t)  we  obtained  in  Eq.  (3.9),  particularly  removing  the  drift  and  obtaining  the  process 
in  Eq.  (3.10),  where  only  the  version  of  the  Brownian  motion  changes. 

This  intuitively  appealing  interpretation  of  a  particle/pixel  motion  in  the  process  of  a  diffusion 
is  shown  in  the  next  section  to  be  particularly  useful  and  insightful  for  developing  more  general  and 
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Figure  3.4:  Middle  row:  Generator  of  symmetrie  random  walk  on  T  is  shown  to  produee  similar 
results  with  those  in  Bottom  Row:  Geometrie  Heat  flow.  The  speeds  of  the  two  algorithms  are 
different.  The  level  set  funetion  is  on  a  191  x  221  grid,  5t  =  0.25. 

feature/shape  adapted  flows. 

3.4  A  New  Class  of  Flows 

The  insight  gained  from  the  tangential  Brownian  motion  on  a  eurve  together  with  the  normal  angle 
0{t,x),  leads  to  the  idea  of  eonstraining  the  Brownian  motion  at  some  speeifie  orientation  angles 
at  eaeh  point  x  .  A  natural  modifieation  of  the  geometrie  heat  equation,  based  upon  the  stoehastie 
framework  presented  in  Seetion  3.3,  is  to  eonstruet  an  SDE  weighted  by  a  earefully  ehosen  fune- 
tional  h{9'^),  {h{.)  G  (^“(M”))  designed  to  eapture  speeifie  features  in  an  image,  and  we  write 
loeally 


dx  (t)  =T^{x  (t))  (c^(x  (t))h(e^{x  (t)))  dt  +  V2  h(e^{x  (t)))  dB(t)). 

Here,  again,  negleeting  the  drift  motion  and  eoneentrating  on  pure  diffusion,  the  Brownian  motion 
in  the  tangent  direetion  is  being  further  eonstrained  at  some  speeifie  orientation  values,  i.e.  at  the 
zeros  of  the  h{0^)  funetion, 

dX{t)  ^V2T^{X{t))h{e^{X{t)))dB{t)).  (3.12) 

Constraining  the  diffusion  of  partieles  at  points  with  speeified  orientations  is  aimed  at  extraeting  de¬ 
sired  features  of  a  eontour  as  it  is  being  smoothed.  Sueh  models  are  generated  by  the  following  elass 
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of  PDE’s,  which  directionally  modify  the  geometric  heat  flow  (3.2),  and  in  this  sense,  generalize  it 
by  making  the  local  generator  of  the  diffusion  SDE  (3.12)  conceivably  arbitrarily  selective. 


du{t,  X ) 

di 


(3.13) 


When  applied  to  an  image,  this  flow  induces  the  following  curve  evolution  on  each  iso-intensity 
contour  C 


dC 


h^{e)KN . 


(3.14) 


3.4.1  Well-posedness  of  the  generalized  model 

Proposition  2  The  corresponding  PDE’s  (3.13)  are  well-posed. 


Proof  2  The  geometric  heat  equation  which  corresponds  to  the  simplest  case  of  this  class  with 
h?{6{t,x  ))  =  1,  Vi,  'ix,  has  been  shown  to  be  well-posed,  and  its  existence  and  uniqueness 
properties  may  be  found  in  [6,  31,  38].  The  operator  of  the  geometric  heat  equation  is  given  by 

f)qj 

L[u]  =  C[u]  -  —  =  0  (3.15) 


where 


£[«]  =  E 


d‘^u 

dxi  dxj 


sin 


^  0  Uxx  —  2  sin  9  cos  9  Uxy  cos^  9  u 


yyi 


(3.16) 


is  the  principal  part  of  the  operator  L.  The  matrix  of  coefficients  [aif  is  positive  semi-definite  with 
the  eigen  values  1  and  0.  If  we  multiply  this  matrix  by  a  positive  function,  it  remains  positive  semi- 
definite.  Such  elliptic-parabolic  operators  satisfy  a  maximum  principle  (see,  for  example,  [119]). 
In  our  case  we  multiply  by  a  non-negative  function  h?{9)  which  can  be  made  strictly  positive  by 
adding  a  very  small  number,  e  >  0, 


[h‘^{9)  e]Cu  >  0. 

This  results  in  a  family  of  nonlinear  parabolic  equations  each  of  which  satisfies  a  strong  maximum 
principle.  Our  operator  is  obtained  in  the  limit  as  e  — )■  0. 
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3.4.2  Polygon  yielding  diffusions 


The  geometric  heat  equation  is  a  rotationally  invariant  flow  which  evolves,  as  mentioned  earlier,  any 
shape  into  a  circle  [54].  It  is  the  only  rotationally  invariant  shape  evolution  in  Euclidean  space.  If 
we  wish  to  capture  more  general  shapes  (triangles,  squares,  etc. . . )  it  is  only  then  natural  to  consider 
flows  which  are  not  rotationally  invariant.  Such  a  class  is  given  by  the  form  (3.14)  when  h{0)  is 
chosen  to  be  other  than  a  constant.  If  we  are  particularly  interested  in  polygons,  we  may  consider 
periodic  functions  (whose  periodicity  is  dictated  by  the  desired  number  of  vertices)  such  as 


HO) 


cos {n6) 
sin  {n6), 


leading  to  curve  evolution  equations  of  the  form 


dC 


cos^  {n9)KN 


or 


dC 


sin^  {n9)KN . 


(3.17) 


(3.18) 


If  we  apply  (3.18)  to  a  convex  shape,  there  will  be  2n  points  on  the  curve  which  do  not  diffuse 
(corresponding  to  the  zeros  of  cos{n9)  or  sin(n0))  at  equally  separated  rotations  of  the  unit  normal 
N .  As  the  unit  normal  moves  further  and  further  away  from  these  angles,  the  diffusion  increases. 
It  hence  makes  sense  to  expect  a  curve  to  develop  vertices  (points  of  maximal  curvature)  at  these 
points. 

Lemma  1  The  angle  of  a  unit  normal  does  not  change  at  points  where  the  chosen  function  hf{9) 
vanishes.  Those  points,  in  turn,  remain  fixed  for  a  short-time,  and  their  speed  remains  at  zero. 

Proof  3  Assume  that  a  family  of  curves  C  (t,  p),  where  p  is  any  parameter  along  the  curve,  evolves 
according  to  the  evolution  equation 


dC 


a{t,p)T  -\-P{t,p)N 


(3.19) 


The  evolution  equation  for  the  angle  of  the  unit  normal  is  given  in  [76]  as 

89  —  1 

—  =  —  [/3p-  aKg]  (3.20) 

dt  g  ^ 

where  g  =  \\C p\\  =  is  the  length  along  the  curve  (metric).  If  we  consider  the  case 

a  =  0  and  /3  =  —H{9)k  (following  the  convention  used  by  the  authors  in  [76]),  which  corresponds 
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to  the  form  of  the  deformation  we  proposed,  the  orientation  evolution  is  governed  by 

=  ^{2h{e){h{e))pn  +  h^{e)np} 

=  ^h{e)  {2{h{e))p  K  +  h{e)  Kp]  (3.21) 

Notice  that  ^  =  0/or  those  points  at  which  h{9)  =  0. 

We  note  that  in  [76],  the  orientation  of  a  curve  is  defined  as  the  angle  subtended  by  the  tangent 
and  the  x-axis,  whereas  here  we  define  0  as  fhe  angle  subfended  by  fhe  normal  and  fhe  x-axis.  There 
is,  however,  a  complefe  equivalence  in  so  far  as  fhe  evolufion  equation  of  fhe  angle  0  is  concerned. 

In  lighf  of  fhe  above  developmenf,  we  can  fhus  sfafe  fhaf  fhe  zeros  of  fhe  function  h{0)  lead  fo 
fixed  end  poinfs  of  curve  segmenfs.  Fixing  fwo  end  poinfs,  say  ai  and  02  ,  we  examine  fhe  evolufion 
of  curvafure,  whose  general  form  is  given  by  (in  [76]) 

Ok  d‘^13  Ok  2 

where  s  is  fhe  arc-lengfh  paramefer  along  fhe  curve.  When  subsfifufing  a  =  0  and  f3  =  —hf{9)K 
info  Ibis  equafion,  we  have 

^  =  [h\9)KU  +  h\9)K^ 
c)k 

—  =  [(hH9))ssi^  +  2(hH9))sns+hH0)i^ss]  +  hH0W 

c)k 

—  =  1^(0^  +hH0)  +  {hH0))sj  ^  +  HhH0))s 

diffusion  ferm  reaction  ferm 

This  clearly  demonsfrafes  fhaf  a  regularizing  diffusion  lakes  place,  since  fhe  mulfiplicafive  faclor 
h?{0)  never  becomes  negafive  (which  would  resulf  in  an  ill-posed  backward  diffusion).  In  addifion, 
we  have  fhe  reaction  term  which  is  composed  of  funclions  of  k,  nf,  and  Kg. 

We  have  hence  shown  fhaf  wifh  fixed  end  poinfs,  a  particular  curve  segmenf  subjecf  fo  fhe  new 
evolution  equation  for  fhe  curvafure  shown  above,  resulls  in  a  slraighl  line  as  a  final  solution. 

Now,  we  can  slate  a  Iheorem  where  we  pul  our  argumenl  of  convergence  fo  regular  polygons. 

Theorem  1  A  convex  curve  C  subject  to  the  evolution  Ct  =  h?{0)KN  will  converge  to  an  M- 
sided,  regular  polygon  whose  M  vertices  will  be  formed  at  those  vanishing  points  of  the  function 
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h^{e). 

The  proof  of  this  theorem  ean  be  eompleted  using  the  are-length  evolution  equation 

where  ds  denotes  the  incremental  arclength  of  C .  Since  the  integrand  is  strictly  positive,  we  see 
that  a  curve  will  continue  to  shrink  until  curvature  vanishes,  that  is  the  curve  segment  converges 
to  a  straight  line  between  the  end  points  ai  and  02 ■  This,  in  conjunction  with  the  above  lemma 
completes  the  proof  of  the  theorem.  ^ 

3.5  Experimental  Results 

3.5.1  Examples  in  Polygonization 

To  substantiate  the  stated  theorem,  and  to  intuitively  illustrate  our  flows,  we  next  present  simulation 
results.  In  our  experiments  with  contours,  we  use  the  narrowband  implementation  of  the  level 
set  method  developed  in  [2].  The  time  step  is,  5t  =  0.2.  Starting  with  a  circular  shape,  the  flow 
Ct  =  h?{9)KN  evolves  it  towards  a  specific  polygon,  i.e.  it  produces  an  n  —  pone  shape  depending 
on  the  specific  function  h{9).  Several  examples  on  morphing  of  a  circle  into  different  shapes  are 
shown  in  Fig.  3.5.  This  is  one  potential  application  of  the  proposed  flows  in  computer  graphics. 


h{9)  =  sin(1.50)  h{9)  =  sin(20)  h{9)  =  sin(2.50)  h{9)  =  sin(40) 

Figure  3.5:  Morphing  of  a  circle  into  different  shapes  by  the  given  flows  is  demonstrated. 

where  the  ability  to  morph  a  shape  into  a  known  other  shape  with  an  efficient  algorithm  is  required 
for  numerous  applications.  In  addition  to  illustrating  the  propagation  of  the  proposed  flows  in 
several  snapshots.  Fig.  3.5  also  provides  a  quantitative  and  an  objective  means  for  characterizing 

'Note  added  in  Proof:  We  recently  found  out  per  Dr.  Osher  at  UCLA  that  Peng,  Osher,  Merriman  and  Zhao,  [116], 
have  also  independently  proposed  flows  similar  to  those  described  in  this  thesis,  albeit  from  a  totally  different  perspective, 
and  with  a  convective  trend. 
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the  performance  of  these  algorithms  in  preserving  corners.  It  can  be  observed  that  the  position  of 
the  desired  feature  locations,  i.e.  the  orientations  at  which  the  vertices  of  the  final  polygon  are 
to  be  formed,  are  well  preserved.  There  may  inevitably  be  one-to-two  pixel  displacements  due  to 
numerical  implementation  effects,  on  account  of  the  finite  precision  of  the  computations  and  the 
finite  resolution  of  the  grid  (which  affects  almost  all  image  processing  algorithms). 

The  flows  are  also  applied  to  a  variety  of  convex  shapes  shown  in  part  (a)  of  each  figure:  Fig.  3.6 
and  Fig.  3.7.  In  Fig.  3.6,  the  shapes  in  part  (b)  were  obtained  by  using  h?{9)  =  cos^  {29)  via  the 
following  curve  evolution 

Ct  =  cos‘^{29)KN ,  (3.22) 


while  the  shapes  in  part  (c)  were  obtained  using  h?{9)  =  sin^  (20).  In  both  cases,  we  expect  to 
obtain  four-sided  and  regular  polygons.  The  zeros  of  cos (20)  and  the  zeros  of  sin(20)  are  however 
45  degrees  out  of  phase.  As  such,  we  see  the  evolved  shapes  in  part  (b)  taking  the  form  of  a  square, 
whereas  the  evolved  shapes  in  part  (c)  take  the  form  of  a  diamond,  corresponding  to  a  45  degree 
rotation  of  the  shapes  in  part  (b).  In  Fig.  3.7,  we  see  the  effect  of  using  different  periods.  The  shapes 
in  part  (b)  are  obtained  using  C  t  =  cos^(30)KiV ,  while  the  shapes  in  part  (c)  are  obtained  using 
C  t  =  sin^(1.5(0  —  'k/2))kN  .  In  the  first  case,  we  expect  6  vertices,  and  in  the  second  case  we 
expect  3  vertices.  Our  expectations  match  the  results  shown  in  part  (b)  and  (c),  where  we  observe 
hexagonal  and  triangular  shapes,  respectively. 


(c) 


Figure  3.6:  (a)  Initial  set  of  shapes  (b)  Flow  C t  =  cos‘^{29)kN  (c)  Flow  Ct  =  sin^(29)KJV  . 


Figure  3.7:  (a)  Initial  set  of  shapes  (b)Flow  C  t  =  cos"^  {W)kN  ,  which  tends  to  produce  hexagons, 
(c)Flow  C t  =  sin^(1.5(0  —  'k/2))kN  ,  which  tends  to  produce  triangle-like  shapes. 


These  two  figures  also  suggest  an  important  potential  application  of  the  proposed  flows,  namely 
shape  recognition.  A  typical  scenario,  would  consist  of  an  unknown  shape,  which  is  evolved  with 
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one  of  the  polygonizing  flows  whose  parameters  are  known  a  priori.  The  evolution  whieh  subse¬ 
quently  results  in  the  least  ehange  in  the  input  shape  reveals  the  elosest  shape  eategory  the  test  shape 
may  belong  to. 

Reeall  that  we  may  also  apply  these  flows  to  the  level  sets  of  an  image  in  the  same  manner  that 
the  geometrie  heat  equation  may  be  applied.  This  gives  rise  to  a  family  of  anisotropie  smoothing 
filters  whieh,  unlike  the  geometrie  heat  equation,  are  not  rotationally  invariant.  This  feature  ean 
be  useful  in  smoothing  noisy  images  where  eorners  and  edges  are  priorly  known  to  have  eertain 
orientations.  These  diffusions  are  effeeted  by  PDE’s  of  the  following  form: 


ut  =  h^{e)v  ■ 


(3.23) 


Note  that  the  trigonometrie  expressions  we  have  eonsidered  for  h?{9)  ean  be  written  in  terms  of  the 
first  derivatives  of  u,  for  example 


cos^(20) 


[ul  +  ul) 


2\2 


and  sin^(20) 


(2UxUy) 


allowing  one  to  implement  the  PDE  without  having  to  eompute  the  orientation  of  the  unit  normal 
to  eaeh  level  eurve.  Note  that  these  expressions  involve  only  first  order  derivatives  and  therefore  do 
not  alter  the  quasi-linear  strueture  of  these  seeond  order  flows. 

The  intended  applieation  of  the  proposed  flows  in  this  ehapter,  i.e.  the  smoothing  of  struetures 
along  the  orientation  of  salient  lines  in  both  eurves  and  images  will  be  illustrated  in  the  next  two 
subseetions,  respeetively. 


3.5.2  Examples  in  Feature-Preservation 

Eeature-preserving  properties  as  well  as  polygonal  approximation  properties  of  the  proposed  flows 
will  be  demonstrated  in  this  seetion.  We  illustrate  the  idea  of  eapturing  different  polygonal  features 
of  shapes  by  our  proposed  flows  on  the  following  examples. 

The  first  example  is  a  “ehef”  shape  with  both  round  and  polygonal  features  as  shown  in  Eig.  3.8. 
The  geometrie  heat  flow  C  t  =  kN  ,  {n  =  0),  evolves  these  features  into  eireles  as  shown  in  the 
seeond  row  of  Eig.  3.8  for  time  points  t  =  40,80, 160.  Partieularly,  at  f  =  160,  most  parts  of 
the  shape  turns  into  ineomprehensible  blob-like  struetures.  In  eontrast  to  this,  polygonal  features 
of  the  ehef  like  his  nose,  and  tray,  are  preserved  by  the  flow  C  t  =  sm‘^{20)KN ,  {n  =  2),  whieh 
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favors  diamond-like  structures  (see  third  row  of  Fig.  3.8  for  t  =  40,  80, 160).  Similarly,  the  flow 
C  t  =  cos"^  {46) kN  ,  {n  =  4),  favors  octagonal  features  as  shown  on  the  fourth  row  of  Fig.  3.8, 
which  is  observed  at  chef’s  hat  at  all  time  points  t  =  40,  80, 160.  The  regularity  of  these  flows  is 
readily  observed  through  the  smoothness  of  the  resulting  shapes.  When  we  view  each  row  from  left 
to  right,  we  observe  a  progression  from  finer  to  coarser  scale.  The  scale-spaces  produced  by  our 
modified  flows  in  the  last  two  rows  are  visually  more  pleasing  since  corners  are  preserved,  whereas 
in  the  row  above  we  see  them  smoothed  away  by  the  pure  geometric  heat  flow. 


Actual  Shape 


t=80 


Figure  3.8:  Each  row  corresponds  to  a  curve  evolution  method  with  different  n,  1**  row:  C  t  =  kN  , 
2”*^  row:  C  t  =  sin^  {29)kN  ,  3’'*^  row:  C  t  =  cos^  {4:d)KN . 


The  second  shape  example  is  a  fish  which  contains  some  fine  detail  structures  as  well  as  coarse 
features  (Fig  3.9).  The  second  row  shows  the  result  of  the  geometric  heat  flow  C  t  =  kN  ,(n  =  0), 
which  smoothes  away  not  only  fine  features  but  some  coarse  features  as  well  (the  fins  for  example). 
The  results  of  the  flow,  C  t  =  cos‘^{26)kN  ,{n  =  2),  are  shown  in  the  third  row  of  Fig  3.9.  In 
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this  case,  rectangular  features  are  preserved  for  longer  periods  throughout  the  evolution.  Finally, 
the  flow,  C  t  =  siTo?{4:6)KN ,  {n  =  4),  is  depicted  in  the  last  row,  preserving  octagonal  features  as 
shown  in  the  nose  and  the  dorsal  fins. 


Actual  Shape 


Figure  3.9:  Each  row  corresponds  to  a  curve  evolution  method  with  different  n,  1**  row:  C  t  =  kN  , 
2”*^  row:  C  t  =  cos^  {29)kN  ,  3’’'^  row:  C  t  =  sin^  (40)kA/’  . 


In  the  third  example,  we  start  with  a  noisy  shape  at  time  t  =  0  shown  in  Fig.  3.10.  This  shape  is 
evolved  with  the  geometric  heat  flow  C  t  =  kN  ,  (n  =  0),  the  flow  C  t  =  sin^(1.5(0  —  'k/2))kN  , 
(n  =  1.5),  the  flow  Ct  =  s\v?{29)nN ,  {n  =  2),  and  the  flow  Ct  =  cos^(2.5(0  —  n/2))KN , 
(n  =  2.5),  as  shown  in  Fig.  3.10.  The  geometric  heat  flow  at  the  top  row  quickly  smoothes  corners 
of  the  shape  out,  and  at  coarser  scales,  the  shape  loses  all  of  its  features.  The  initial  shape  converges 
to  a  circle  in  spite  of  the  global  feature  of  the  plane  being  a  polygonal  shape.  This  motivates  the 
application  of  the  geometric  heat  flow  with  a  sin^(n0)  factor,  where  n=1.5,  and  whose  weak  limiting 
shape  is  a  triangle  which  intuitively  matches  the  coarser  form  of  the  given  plane  shape.  Similarly, 
forn  =  2  and  n  =  2.5,  different  features  of  the  shape  are  preserved,  and  persist  over  a  much  longer 
time  period  as  can  be  observed  from  the  column  of  shapes  at  f  =  400.  Note  that  the  geometric  heat 
flow  result  at  the  top  quickly  washes  out  any  similarity  to  the  actual  shape,  whereas  the  results  of 
the  other  three  flows  preserve  the  global  shape  as  well  as  some  finer  details  on  the  wings,  the  tail, 
and  the  head  part. 
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3.5.3  Examples  with  Grayscale  Images 

The  proposed  flows  may  also  be  applied  to  images  in  a  straightforward  fashion.  For  the  ease 
h?{0)  =  cos^(20),  all  level  sets  of  the  image  are  driven  to  reetangles,  thereby  enhaneing  those 
features  in  an  image.  Sueh  features  ean  be  found  in  eontemporary  buildings  where  one  example  in 
NCSU,  Centennial  Campus  is  shown  in  Fig.3. 11(a).  The  part  of  the  building  image  with  an  addi¬ 
tive  Gaussian  noise  is  shown  in  Fig.3. 11(b).  The  2”*^  row  shows  the  results  of  the  geometrie  heat 
flow  ut  =  at  t  =  10, 20,  50.  The  noisy  image  at  f  =  0  is  smoothed  out  very  quiekly  at  the 
expense  of  rounding  off  all  the  eorners  beeause  the  level  sets  of  the  image  eonverge  to  eireles.  The 
yd,  j.Q^  shows  the  ut  =  cos^(20)u^^  flow  results  at  the  same  time  points  t  =  10,20,50.  Sinee 
the  diffusion  is  eonstrained  in  order  to  drive  image  level  sets  to  reetangles,  the  removal  of  noise  is 
slower.  However,  the  reetangular  features  still  nieely  appear  after  noise  removal  (see  the  image  on 
the  right),  making  it  worthwhile  to  slow  down  the  denoising  effeet  of  the  geometrie  heat  flow  as 
deemed  appropriate. 

In  Fig.  3.12,  an  experiment  involving  diamond-like  shapes  in  the  image  taken  from  a  poster  on  a 
wall  is  shown.  In  the  middle  row,  rounding  effeets  on  diamond  shapes  performed  by  the  geometrie 
heat  flow  are  elearly  observed  during  the  evolution  of  this  image.  The  proposed  flow,  shown  in 
the  bottom  row,  takes  the  form  ut  =  sin^(20)u^^  for  this  partieular  shape,  and  partieularly  adapted 
to  earrying  out  a  shape-based  smoothing  whieh  takes  plaee  at  the  boundaries  of  the  diamonds. 
The  slight  blurring  effeet  on  the  pieture  at  eontinued  applieation  however  is  due  to  the  interaetion 
between  eonseeutive  level  eurves. 

A  photograph  taken  by  pathfinder  in  mars,  shown  in  Fig.  3.13,  is  argued  to  be  a  hexagon-shaped 
structure  on  mars’  surface.  The  particular  flow  adapted  to  this  shape  is  given  by  ut  =  cos^(30)u^^, 
and  the  resulting  images  at  the  second  column  of  the  figure  demonstrate  a  better  smoothing  perfor¬ 
mance  at  the  boundaries  of  the  hexagon  when  compared  to  the  images  in  the  first  column  processed 
by  the  geometric  heat  equation.  From  low  scales  to  very  large  scales,  the  hexagon-adapted  flow 
enhances  and  keeps  on  highlighting  the  related  structure. 

A  noise  contaminated  Aerial  image  is  shown  in  Fig.  3.14(a).  The  geometric  heat  equation  (see 
2^d  ^  _  20, 40,  80)  sweeps  away  the  shape  information  of  the  important  details  such  as  the  city 

on  the  left  bottom,  the  white  bright  rectangle  on  the  right  bottom,  and  the  black  feature  at  the  top. 
The  three  images  resulting  from  the  ut  =  cos^(20)  flow,  are  quite  sharp  at  the  edges  between 
both  low  and  high  contrast  fields,  therefore  more  useful  in  recognition  of  details  as  well  as  removing 
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noise. 


A  last  example  is  shown  in  Fig.  3.15,  where  windows  and  a  roof  of  a  seetion  of  a  house  are  seen. 
On  the  left  eolumn,  the  results  of  the  geometrie  heat  flow  ut  =  at  times  t  =  40,  80, 160,  are 
shown,  and  on  the  right  those  of  ut  =  cos^(20)i/^^  at  similar  time  steps.  The  noise  is  sueeessfully 
removed  by  the  geometric  heat  equation  whose  smearing  effect  on  different  regions  into  one  another 
is  also  slow,  at  a  cost  of  a  problematic  rounding  off  of  corners.  At  time  t  =  160  for  the  result  on  the 
right  bottom,  approximately  the  same  amount  of  noise  as  that  of  geometric  heat  equation  at  f  =  40 
is  removed,  and  in  addition  to  that  the  corners  are  still  well-preserved. 

3.6  Conclusions 

In  this  chapter,  also  in  [141],  we  have  formulated  a  local  stochastic  view  of  a  nonlinear  filtering 
technique,  namely  the  geometric  heat  equation.  The  theory  of  stochastic  differential  equations  pro¬ 
vides  a  microscopic  view  of  a  system,  and  through  a  local  linearization  of  the  nonlinear  geometric 
heat  equation  we  have  provided  an  alternative  macroscopic  view  of  this  equation.  We  then  mod¬ 
ified  this  macroscopic  description  to  propose  new  flows  that  vanish  at  pre-defined  directions.  We 
showed  that  these  flows,  although  rotationally  non-invariant,  are  capable  of  smoothing  along  priorly 
known  orientations  of  salient  lines  in  both  curves  and  images,  leading  to  preservation  of  polygonal 
structures.  In  the  context  of  curve  evolutions,  curves  evolved  with  the  new  flows  are  morphed  into 
a  limiting  polygonal  shape,  this  approach  may  hence  be  relevant  to  shape-morphing  applications  in 
computer  graphics.  Another  application  of  these  polygonizing  flows  is  in  classification  of  a  shape 
after  its  filtering  via  a  certain  set  of  these  flows,  and  its  identification  according  to  the  outputs  of 
these  different  filters. 

Our  main  idea  of  stopping  the  diffusion  with  a  function  of  the  normal  angle  0,  has  an  important 
characteristic  that  h{-)  function  should  be  periodic  since  the  angle  0  lives  on  the  unit  circle.  We  have 
presented  in  [140, 142],  a  variant  of  this  idea  where  the  modification  of  the  geometric  heat  equation 
is  accounted  by  a  function  that  depends  on  the  unit  normal  angle  difference  between  neighboring 
points  on  the  tangent  line  at  a  point  of  the  curve.  This  also  leads  to  a  feature-driven  smoothing 
of  curves,  (with  different  characteristics  than  the  flows  presented  in  this  chapter),  which  has  been 
used  in  filtering  level  sets  of  images,  particularly  of  synthetic  aperture  radar  imagery,  and  also  in  the 
target  recognition  of  the  latter. 
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t=2 


t=5 


t=10 


t=20 


Figure  3.12:  (Top)  Diamonds  image  (Middle  row)  Geometrie  heat  flow  ut  =  (Bottom  row) 
Flow  Ut  =  sin^(20)  u^^. 
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Figure  3.13:  (Top)  An  image  from  mars  pathfinder,  (First  eolumn)  Geometrie  heat  flow  ut  = 
(Seeond  eolumn)  Flow  ut  =  cos^(30)  (Image:  Origin  NASA,  exposed  by  and  eourtesy  of  N. 
Coombs  &  NPAAG  1998.) 
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Figure  3.14:  Top:  Aerial  image;  2”*^  row:  geometrie  heat  flow  ut  =  (left  to  right)f  =  20, 40,  80; 
yd'  row:  flow  ut  =  cos^  {20)u^^,  (left  to  right)  t  =  20, 40,  80. 
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Figure  3.15:  Top:  House;  Left:  Geometrie  heat  flow  ut  =  (top  to  bottom)  t  =  40,80, 160; 
Right:  Flow  ut  =  cos^(20)  (top  to  bottom)  t  =  40,  80, 160. 
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Chapter  4 


Unsupervised  Texture  Segmentation  by 
Information-Theoretic  Active  Polygons 


This  chapter  introduces  new  ideas  for  image  segmentation  problem  through  new  polygon-propagating 
equations  in  the  active  contour  framework.  Furthermore,  adoption  of  an  information-theoretic  cri¬ 
terion  as  an  energy  functional  of  an  active  polygon  is  described  for  an  unsupervised  texture  seg¬ 
mentation  purpose.  An  overview  on  the  general  segmentation  problem  was  given  in  Chapter  2.3. 
After  briefly  reviewing  this  literature  in  Section  4.1,  an  overview  on  texture  analysis  and  modeling 
is  presented  in  Section  4.2,  along  with  our  motivations,  strategy,  and  the  advantages  of  pursuing  an 
active  polygon  model  over  adapting  an  active  contour  model.  In  Section  4.3,  we  derive  the  ODEs 
to  obtain  motion  equations  for  polygon  vertices.  In  Section  4.4,  we  introduce  Jensen-Shannon  cri¬ 
terion  for  evolution  of  a  single  contour  and  multiple  contours  with  the  goal  of  texture  segmentation. 
In  Section  4.5,  a  novel  polygon  regularizer  to  avoid  degeneracy  during  a  polygon  propagation  is 
introduced.  To  our  knowledge,  this  type  of  global  geometric  regularizer  has  not  been  used  even 
in  related  spline-based  active  contour  models.  Results,  conclusions  and  discussions  are  given  in 
Section  4.6. 


4.1  Introduction 

The  problem  of  image  segmentation  refers  to  the  partitioning  of  a  given  image  domain  into  regions 
which  are  distinct  in  some  sense.  These,  in  turn,  are  expected  to  correspond  to  meaningful  parts  of 
objects  in  the  image.  Image  segmentation  is  generally  viewed  as  an  essential  first  step  in  low  level 
vision  and  as  providing  a  mechanism  for  an  automatic  analysis  of  image  contents.  Some  important 
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applications  include  automatic  target  recognition,  remote  sensing,  robot  vision  and  guidance,  auto¬ 
matic  visual  inspection  in  manufacturing  processes,  biomedical  image  analysis,  tracking  objects  in 
motion.  Segmentation  is  generally  based  on  first  determining  features  which  delineate  and  reliably 
distinguish  different  regions.  Some  common  features  which  discriminate  the  region  characteristics 
are  intensity,  color,  and  texture.  Selection  of  a  “good”  criterion  is  critical  to  capturing  and  evaluating 
the  features  which  yield  a  partition  of  the  image  domain  into  distinct  regions. 

Earlier  image  segmentation  methods  were  based  on  thresholding  [1 12, 134]  and  on  local  filtering 
techniques  [22, 96, 97],  which  make  use  of  local  information,  and  rely  on  high  image  gradient  values 
to  detect  prominent  boundaries  between  regions.  This  in  turn  makes  them  sensitive  to  noise  and 
affects  the  continuity  of  the  edge  contours  as  a  result.  Region  growing  techniques  [100]  partition  an 
image  domain  by  merging  regions  according  to  some  statistical  criterion  such  as  region  variances, 
but  their  disadvantage  is  that  they  often  generate  irregular  boundaries  [153].  Snake  [70]  and  active 
contour  methods  define  energy  funcfionals  whose  local  minima  comprise  a  sef  of  solutions,  e.g. 
boundaries  of  regions.  The  fwo  main  sfreams  of  fhoughf  are  as  follows: 

•  Inspired  by  fhe  celebrafed  geomefric  heaf  equation  [54],  which  may  be  obfained  as  a  gradienf 
descenf  flow  of  a  Euclidean  arc  lengfh,  geomefric  acfive  confour  models  [23, 25, 72, 93, 130], 
were  developed  as  gradienf  flows  of  a  modified  Euclidean  arc  lengfh.  A  modifying  weighf 
which  is  almosf  zero  near  edges,  and  almosf  one  when  if  is  far  from  fhe  edges,  can  be  chosen 
fo  be  inversely  proportional  fo  fhe  magnifude  of  fhe  image  gradienf.  These  models  are  hence 
boundary-based,  and  are  only  sensitive  fo  dafa  near  fhe  curve  (very  local).  This  makes  fhem 
prone  fo  noise  variability  and  fo  inifial  confour  placemenf. 

•  To  overcome  fhis  problem,  region-based  active  confours,  which  use  bofh  local  and  global 
information  were  proposed  [26, 27, 29, 121, 123, 149, 153].  These  models  mainly  assume  fhaf 
fhe  image  consisfs  of  a  finite  number  of  regions,  fhaf  are  characferized  by  a  pre-defermined 
sef  of  fealures  or  sfafislics  such  as  means,  and  variances.  An  energy  funclional  is  consfrucfed 
fo  pull  fhese  sfafislics  aparl.  One  advanlage  over  fhe  geomefric  active  confours  is  fhaf  Ihere  is 
no  need  fo  calculate  image  gradienfs  which  are  usually  sensitive  fo  noise,  albeif  af  a  cosf  of 
additionally  imposed  assumpfions  on  fhe  image. 

The  popularity  of  fhese  active  confour  models  was  parficularly  due  fo  an  efficienl  implemenfafion 
scheme  via  partial  differenlial  equations  (PDEs),  referred  fo  as  level  sef  mefhods  developed  by  Osher 
and  Sefhian  for  propagating  interfaces  [111,  127]. 
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Global  optimization  approaches  to  image  segmentation  are  based  on  energy  functions  minimiz¬ 
ing  different  criteria  such  as  Minimum  Description  Length  (MDL)  [87],  and  Bayesian  [17].  Other 
methods  were  more  hybrid  in  nature,  such  as  the  Region  Competition  method  of  Zhu  and  Yuille 
in  [153],  which  adopt  a  generalized  MDL  criterion  to  essentially  include  the  costs  for  coding  the 
intensity  of  each  pixel  inside  a  region  according  to  a  prespecified  probability  distribution  whose 
parameters  are  to  be  estimated.  A  region  merging  step  is  subsequently  carried  out  while  ensuring 
a  decreasing  total  energy.  A  well-known  mathematical  model  proposed  by  Mumford-Shah  [103] 
combined  a  boundary  extraction  goal  via  an  assumption  of  approximating  an  image  by  smooth 
functions  in  each  region.  Level-sets  active  contour  implementation  of  this  model  have  recently  been 
introduced  in  [28, 137].  In  fact  region-based  active  contour  models  mentioned  in  the  previous  group 
can  also  be  considered  in  this  group,  because  their  goal  is  also  to  extract  boundaries  of  homogeneous 
regions  within  an  image  by  minimizing  an  energy  functional. 

A  particularly  powerful  cue  in  visual  perception  is  the  set  of  textural  features.  Texture  segmen¬ 
tation,  which  is  the  task  of  parsing  the  image  domain  into  a  number  of  regions  such  that  each  region 
has  the  same  textural  properties  is  a  challenging  problem  [68,94].  The  problem  of  unsupervised 
texture  segmentation,  which  proceeds  without  a  priori  information  about  the  textural  characteristics 
of  objects  in  a  given  image,  is  particularly  challenging  and  remains  largely  an  open  research  issue 
in  computer  vision.  This  along  with  its  ubiquitous  applications  motivate  further  investigations  into 
developing  refined  and  sophisficafed  mefhods  fo  accurafely  reflecf  and  discriminafe  fexfured  regions 
wifhouf  supervision. 

For  fexfured  images,  defining  a  homogeneify  measure  is  a  mafhemafically  challenging  prob¬ 
lem.  We  are  going  fo  propose  a  new  fexfure  segmenfafion  approach  which  is  an  efficienf  improve- 
menf  over  fhe  existing  class  of  mefhods.  Specifically  we  propose  an  informafion-fheorefic  measure, 
namely  Jensen-Shannon  divergence,  which  fo  our  knowledge  has  nof  been  used  in  fhe  active  contour 
framework,  and  seeks  fo  disfinguish  and  separafe  probabilify  densify  funcfions  of  various  regions 
in  an  image.  This  measure  makes  use  of  higher  order  sfafisfics  beyond  mean  and  variance  by  a 
compufafionally  efficienf  manner  which  does  nof  require  histograms. 
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4.2  Motivation  and  Previous  Work 


4.2.1  Related  Work 

For  texture  analysis  and  modeling,  mainly  two  groups  of  approaehes  and  their  eombinations  have 
been  developed: 

•  Visual  system’s  deeomposition  of  retinal  images  into  a  set  of  subbands,  has  been  the  souree 
of  inspiration  for  filtering  teehniques  based  on  filterbanks  of  for  instanee  Gaussians,  and 
Gabor  funetions  to  eharaeterize  different  attributes  of  a  texture  at  different  orientations  and 
seales  [68, 131].  Although  multichannel  filtering  techniques  have  been  shown  to  be  efficient 
at  capturing  local  spatial  features,  problems  such  as  optimal  choices  of  filters,  and  fusion  of 
their  outputs  remain  open. 

•  Statistical  modeling  approaches  characterize  texture  images  as  arising  from  probability  dis¬ 
tributions  on  Markov  Random  Fields  (MRF’s)  [13,39,51,94],  and  hence  mainly  rely  on  a 
Markov  property  where  each  pixel  is  statistically  dependent  only  on  a  certain  set  of  neighbor¬ 
ing  pixels.  The  segmentation  is  achieved  through  a  minimization  of  a  maximum  a  posteriori 
(MAP)  criterion  of  the  observed  image.  An  advantage  of  these  approaches  is  that  they  yield 
the  parameters  of  the  underlying  probability  distribution  which  in  turn  affords  one  an  ability 
to  synthesize  texture  images  by  sampling.  Limitations  of  commonly  used  MRF  models,  on 
the  other  hand,  are  due  to  the  fact  that  only  the  first  and  second  order  statistics  may  tractably 
be  used,  while,  it  is  widely  recognized  that  many  textures  are  strongly  non-Gaussian  regard¬ 
less  of  the  neighborhood  size  [154].  An  analytical  probability  density  for  modeling  clutter  in 
natural  images  was  proposed  in  [55]. 

•  Methods  which  attempt  to  unify  the  afore-mentioned  approaches  have  recently  received  at¬ 
tention.  Zhu  and  Yuille  [153]  applied  their  model  to  texture  segmentation  where  two  lo¬ 
cal  orientations  of  texture  elements  obtained  through  a  Gaussian  filtering  of  image  gradient 
components  are  used  as  a  multivariate  texture  image  in  their  region  competition  model  with 
underlying  multivariate  Gaussian  distributions. 

Zhu,  Wu,  Mumford  [154],  also  developed  a  probability  model  built  on  the  texture  features 
extracted  by  a  set  of  filters  aimed  at  capturing  the  properties  of  the  texture  at  multiple  scales 
and  orientations.  Their  probability  model  is  a  maximum  entropy  distribution,  which  then 
includes  the  whole  marginal  distribution  in  this  sense,  and  is  hence  more  expressive  than 
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classical  MRF’s.  This  enables  them  to  model  non-Gaussian  textures.  We  should  however 
note  that  one  of  their  main  goals  is  texture  synthesis  which  is  a  dijferent  goal  than  ours, 
which  is  a  texture  segmentation.  This  texture  synthesis  goal  also  bears  a  heavy  computational 
burden,  which  we  mean  to  bypass  in  our  proposed  approach. 

Paragios  and  Deriche  [115],  in  a  supervised  texture  segmentation  problem,  also  combined  a 
filtering  approach  with  a  statistical  view  while  accounting  for  boundary-based  (geodesic)  and 
region-based  active  contour  frameworks.  For  a  region-based  attraction  of  the  active  contour, 
a  posterior  probability  of  each  region  with  respect  to  the  texture  filter  output  image  is  used, 
exactly  in  the  same  way  as  in  [153],  with  an  MDL  criterion,  which,  as  well-known  may  also  be 
stated  as  a  MAP  criterion.  Their  boundary  module  is  a  geodesic  active  contour  model  which 
utilizes  an  estimated  boundary  probability  model  to  capture  texture  boundaries  (analogous  to 
capturing  intensity  boundaries  with  high  image  gradients).  In  contrast  to  our  approach,  this 
method  requires  a  prior  knowledge  on  the  existing  textures  in  the  given  image,  and  hence 
a  pre-processing  stage  involving  filter  banks  (e.g.  Gabor),  to  estimate  a  Gaussian  mixture 
probability  distribution  representing  these  pre-defined  textures. 

4.2.2  Synopsis  of  Proposed  Approach 

Our  overall  objective  of  object  (homogeneous  region)  delineation  in  an  observed  image  is  highly 
motivated  by  classification  and  recognition  applications.  This  in  turn  and  in  contrast  to  existing 
active  contours,  motivates  a  parsimonious  and  revealing  representation  such  as  that  based  on  land¬ 
marks/vertices  of  an  approximating  polygon  as  we  propose  in  this  chapter.  In  addition  to  its  parsi¬ 
mony,  this  representation  is  expected  to  accurately  capture  the  prevailing  texture  and  account  for  it 
in  the  course  of  the  analysis. 

Given  an  image,  consisting  of  several  differently  textured  and  simply  connected  regions  (or  for 
our  particular  interest,  textured  objects),  an  approximation  of  their  contours  whose  delineation  is 
sensitive  to  their  textural  characteristics  is  achieved  by  adopting  a  polygonal  model  which  we  de¬ 
velop  in  Section  4.3,  and  by  way  of  an  information-theoretic  measure,  which  we  explore  in  Section 
4.4  for  propagating  polygons.  Combining  these  two  goals  yields  our  information-theoretic -based 
active  polygon  which  is  effective  and  appealing  on  account  of  its  implicit  inclusion  of  higher  order 
statistics  of  the  region  data  and  integrated  along  the  edges. 

To  motivate  our  investigation,  we  begin  with  a  simple  example  of  a  synthetic  image  of  two  dis¬ 
tinct  intensity  regions  where  the  limitation  of  an  active  contour  is  readily  seen  as  a  result  of  excessive 
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Figure  4.1:  Motivational  example  where  a  simple  binary  image  containing  a  quadrilateral  (at  the 
top)  with  additive  uniform  noise  is  to  be  segmented.  Rows2-3-4  show  a  continuous  active  contour 
propagation  with  small,  medium,  large  regularization  respectively.  The  first  two  pictures  in  each 
column  are  the  initial  active  contours  on  the  image,  and  on  a  white  background  separately,  and  the 
last  two  pictures  show  what  can  possibly  be  obtained  best  in  the  continuous  case.  On  the  last  row, 
an  initial  arbitrary  active  polygon  (first  two  columns),  is  propagated  successfully  towards  the  target 
region,  and  the  resulting  polygon  has  only  four  vertices  shown  on  the  fourth  column. 
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additive  uniform  noise,  while  a  vertex-based  polygonal  evolution  as  will  be  deseribed  in  detail  eas¬ 
ily  eaptures  the  underlying  objeet.  This  is  illusttated  in  Fig.4.1.  The  essenee  of  the  aetive  eontour  is 
to  propagate  an  initial  eontinuous  eurve  towards  the  boundaries  between  two  regions  via  a  gradient 
deseent  flow  of  an  energy  funetional  with  both  image-based  external  energy,  and  a  smoothness- 
eonstrained  internal  energy.  Most  aetive  eontour  models,  however,  require  signifieant  regularization 
in  the  presenee  of  signifieant  noise  in  the  image.  With  a  little  regularization  on  the  eontinuous  aetive 
eontour  model,  (first  row),  the  eurve  weaves  around  noise,  with  a  medium  regularization,  (seeond 
row),  still  there  is  no  way  for  an  aetive  eontour  to  get  out  of  loeal  minima  in  this  highly  noisy  ease. 
With  a  very  large  regularization,  the  eurve  starts  to  eonverge  to  the  shape,  but  the  shrinking  effeet 
is  too  powerful,  and  the  eurve  eontinues  to  shrink  without  stieking  to  the  data,  and  eventually  will 
eollapse  to  a  point.  These  three  eases  illustrate  the  eontinuum  of  effeets  from  small  to  large  regu¬ 
larization.  We  should  note  that  this  is  based  on  a  widely  used  arelength-based  regularization,  and 
that  bad  effeets  sueh  as  eorners  being  rounded  off,  and  the  shrinking  effeets  although  not  as  promi¬ 
nent,  will  be  seen  in  a  similar  way  with  rigidity  terms  (penalties  on  squared  eurvature).  The  eontour 
would  again  ignore  the  data,  and  eonverge  to  a  eirele  regardless  of  the  shape  of  the  objeet.  By 
evolving  an  aetive  polygon  with  a  relatively  small  number  of  vertiees,  strong  regularizing  internal 
forees  are  no  longer  neeessary  to  keep  the  eontour  from  “breaking-up”.  A  polygon  laid  on  the  same 
noisy  image  in  Fig.4. 1  (fourth  row),  propagates  towards  the  boundaries  with  a  greater  resilienee  to 
noise,  and  results  in  a  detention  of  the  sharp  eorners.  The  target  objeet  is  well-delineated,  and  the 
resulting  polygon  has  only  four  vertiees,  aeeurately  loeated  at  the  eorners  of  the  objeet.  While  our 
introdueed  model  is  applied  to  a  seemingly  eontrived  simple  example,  nevertheless,  illustrates  the 
propagation  of  an  aetive  polygon  through  noise  free  of  eost  of  sharp  eorners  as  well  as  its  adaptivity 
and  independenee  of  initialization. 

4.2.3  Strategy 

To  highlight  our  development  strategy,  we  deseribe  how  the  objeetives  stated  at  the  outset  may 
be  aehieved  while  elearly  noting  key  eoneeptual  differenees  with  previously  proposed  teehniques. 
The  evolution  of  the  vertiees  of  a  polygon  aeeording  to  an  energy  funetional  is  elearly  eonstrained 
by  their  set  membership  to  the  polygon  and  yields  a  set  of  eoupled  ordinary  differential  equations 
(ODEs)  eaeh  eorresponding  to  a  partieular  vertex.  We  hasten  to  note  here  that  our  approaeh  differs 
in  an  essential  way  from  marker  partiele  methods  [127],  or  earlier  snakes  [70],  in  that  the  number 
of  vertiees  of  our  aetive  polygons  is  mueh  smaller,  and  that  the  motion  of  vertiees  are  obtained 
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by  integration  of  information  over  all  edges  rather  than  motion  of  individual  partieles.  In  marker 
partiele  methods,  vertiees  behave  independently  under  the  influenee  of  image  forees,  and  eoupling 
only  eomes  from  the  internal  forees.  We  note  in  addition  that  with  our  approaeh,  these  vertiees 
naturally  end  up  on  important  features  sueh  as  eorners  on  the  boundary  of  an  objeet.  Rather  than 
marker-partiele  based  teehniques,  our  approaeh  is  philosophieally  mueh  more  related  to,  but  never¬ 
theless  quite  different  from,  spline-based  models,  where  the  entire  eontour  is  represented  by  a  few 
eontrol  points,  sueh  as  [16,45].  To  prevent  degeneraey  during  motion  of  the  vertiees  (e.g.  edges 
interseeting  eaeh  other),  we  propose  a  global  regularizer  teehnique  whieh  uses  eleetrostaties  prinei- 
ples  by  viewing  eaeh  edge  of  the  polygon  as  a  line  of  uniform  eharge.  This  partieular  point  is  also 
of  interest  in  its  own  right  as  it  introduees  a  global  geometrie  regularizer  rather  than  loeal  geometrie 
regularizers  that  have  been  popular. 

To  meet  eaeh  of  the  objeetives  set  out,  we  define  an  information-theoretie  measure  whieh,  based 
on  the  so-ealled  Jensen-Shannon  divergenee  as  an  energy  funetional  unfolds  the  textural  information 
through  the  underlying  probability  distributions  of  the  data,  using  easily  eomputed  image  statisties. 
The  formulation  of  sueh  an  energy  funetional  in  eonjunetion  with  a  polygonal  probing  tool  is  fol¬ 
lowed  by  a  gradient-based  minimization  proeedure  whieh  evolves  the  vertiees  by  way  of  a  eoupled 
set  of  ODEs  just  mentioned. 

Note  that  Jensen-Shannon  divergenee  [90],  with  Shannon  entropy  as  a  speeifie  ease  [36],  is  a 
probabilistie  differenee  measure,  with  properties  of  interest  in  various  applieations.  For  instanee, 
unlike  most  other  divergenee  measures,  it  ean  be  generalized  to  a  finite  number  of  probability  distri¬ 
butions,  whieh  makes  it  suitable  to  an  image  segmentation  seenario  with  several  target  regions.  For 
an  image  registration  applieation,  a  similar  divergenee  measure,  based  on  a  Renyi  entropy  [36],  was 
proposed  by  [56, 58],  referred  to  as  a  Jensen-Renyi  divergenee.  Jensen-Shannon  divergenee  was 
also  applied  to  edge  deteetion  on  regular  or  texture  images  in  [53]  by  measuring  the  probabilistie 
differenees  between  relative  frequeney  of  gray  values  in  two  halves  of  a  small  window  that  is  slid 
over  an  image.  Sinee  it  was  developed  for  edge  deteetion  purposes,  the  result  of  this  method  yields 
uneonneeted  edges,  and  eontours,  and  requires  further  proeessing  sueh  as  eontour  linking  for  a  seg¬ 
mentation  purpose.  It  henee  differs  from  our  approaeh  whieh  aims  at  finding  a  solid  segmenfafion 
of  objeefs  in  an  image  wifhin  a  variafional  framework. 

In  eonfrasf  fo  Zhu  and  Yuille  [153]  for  insfanee,  our  teehnique  does  nol  assume,  for  simplieify 
or  for  any  ofher  purpose,  an  underlying  Gaussian  model  for  fexfure  regions.  Our  approaeh  is  henee 
sfafisfieally  non-paramefrie,  and  image  dafa  drives  fhe  esfimafion  of  fhe  divergenee  befween  fhe 
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texture  distributions  as  will  be  shown.  A  region-snake  in  [32],  is  developed  for  segmentation  of 
images  with  only  two  regions,  a  target  region  and  a  background  region,  which  are  assumed  as 
independent  random  variables,  such  an  assertion,  and  particularly  for  optical  images,  required  a 
prewhitening  on  the  image.  For  segmentation  of  images  of  other  modalities,  such  as  ultrasound 
images  and  SAR  images,  authors  in  [33],  and  in  contrast  to  our  approach,  assume  particular  prior 
densities  which  have  explicit  expressions,  and  then  make  use  of  a  MAP  criterion  minimized  by  an 
iterative  stochastic  optimization  via  a  polygon.  The  technique  of  Paragios  and  Deriche  [115]  was 
developed  for  a  supervised  texture  segmentation,  which  requires  a  full  knowledge  of  existing  texture 
classes/regions,  which  may  be  multiply  connected,  thanks  to  the  level-set  implementation  [111]  of 
the  PDFs.  The  approach  we  propose  differs  and  may  be  viewed  as  a  more  challenging  task,  i.e. 
an  unsupervised  texture  segmentation  scenario,  in  which  the  only  prior  information  is  the  number 
of  different  classes/regions  which  are  assumed  to  be  simply  connected  due  to  our  simple  polygon- 
moving  ODEs. 

Our  approach  in  pursuing  an  active  polygon  model  is  further  advocated  next. 

4.2.4  Advantages  of  Active  Polygons  over  Active  Contours 

Active  polygons  enjoy  some  significant  advantages  and  flexibility  over  existing  techniques: 

•  A  first  advantage  arises  in  texture  segmentation.  The  system  of  ODEs  integrate  local  image 
speed  functions  along  two  adjacent  polygon  edges  when  computing  motion  of  each  vertex. 
Thus,  it  turns  local  image  measurements  at  points  on  the  contour  into  global  information, 
as  these  local  measurements  are  integrated  along  the  polygon  edges.  This  makes  our  algo¬ 
rithm  much  more  reliable  and  robust  in  capturing  texture  boundaries  in  contrast  to  continuous 
curves  which  are  local  in  scope.  Region-based  methods  also  use  global  information  inside 
and  outside  the  curve,  but  their  gradient  flow  incorporates  local  information,  and  pointwise 
on  the  curve,  and  is  hence  not  amenable  to  speed  functions  for  capturing  higher-order  statis¬ 
tics  which  can  not  be  estimated  from  pointwise  measurements,  in  any  case. 

•  A  reliable  detection  of  sharp  object  corners  is  particularly  important  in  recognition  and  classi¬ 
fication  of  man-made  objects,  and  the  accurate  performance  of  our  algorithm  is  very  accurate 
with  no  rounding  effects  at  such  important  feature  points,  sets  it  apart  from  the  continuous 
approaches. 

•  Emerging  multimedia  applications  require  targeted  and  specific  data  in  large  databases,  which 
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in  turn  entails  a  heavy  demand  for  an  effieient  aeeess  and  representation  of  visual  objeets  in 
a  seene.  An  effieient  deseription  of  an  objeet  starts  with  its  eontour  representation,  whieh 
through  a  eontinuous  eurve  requires  a  tremendous  amount  of  information  to  be  stored.  If  a 
eontour  is  represented  by  a  small  number  of  vertiees,  this  representation  may  later  be  eom- 
paetly  stored  and  provide  a  signifieant  and  informative  amount  of  eompression  of  the  multi- 
media  eontent.  There  is  a  vast  amount  of  literature  on  obtaining  polygonal  approximations 
of  eontinuous  eontours  [42,48,  66,  83].  The  goal  is  often  to  reduee  the  amount  of  data  (e.g. 
eomplex  figures  in  earfography,  geographieal  maps),  and  polygons  are  powerful  approxima¬ 
tors  of  shapes  for  use  in  lafer  reeognifion  sfages  [77],  sueh  as  shape  mafehing,  and  shape 
coding.  Several  techniques  for  shape  coding  such  as  [106, 151],  enfailed  firsl  exfracfing  ob- 
jecf  boundaries  and  subsequenfly  finding  a  polygonal  approximation  of  fhe  exlracfed  contour. 
Wifh  goal  of  efficienl  compression,  our  fechnique,  is  hence  advanfageous  in  fhaf  if  capfures 
in  one  shof,  a  polygonal  represenfafion  of  objecfs  in  an  image  avoiding  a  fypically  lengfhy 
search  for  meaningful  locafions,  such  as  high  curvafure  poinfs,  used  in  shape  mafehing.  The 
polygonal  represenfafion  of  an  objeef  from  an  image  is  only  as  good  as  fhe  segmenfafion  of  fhe 
image.  If  is  fhus  more  efficienf  to  invoke  fhe  polygonal  represenfafion  direcfly  in  fhe  course 
of  fhe  image  segmenfafion  process  rafher  fhan  afferwards.  In  shorf,  our  approach  conneefs  fhe 
final  polygonal  represenfafion  we  wish  fo  use  direcfly  to  fhe  image  dafa. 

•  Wifh  a  handful  of  vertices  provided  as  a  resulfing  objeef  description,  our  mefhod  can  ef- 
ficienfly  be  applied  fo  visual  image  refrieval  or  indexing,  and  confenf-based  description  of 
mulfimedia  dafa  which  are  fhe  core  functionalities  of  MPEG-7  sfandard  [30]  . 

•  Anofher  advanfage  of  using  ODEs  rafher  fhan  PDEs  is  a  more  efficienf  implemenfafion.  In 
fhe  discretization,  a  significanlly  larger  lime  step  can  be  chosen  fo  speed  up  fhe  algorilhm, 
whereas  fhe  PDEs  usually  require  very  small  lime  sleps  for  slable  numerical  implemenfafions, 
particularly  when  slrong  regularizers,  which  are  avoided  by  our  model,  are  necessary.  In 
addition,  fhe  identification  of  fhe  interior  and  fhe  exterior  of  a  polygon  is  much  easier  fhan  a 
conlinuous  curve,  Ihereby  providing  additional  speed  up  in  fhe  calculations  of  slafislics  inside 
and  oulside  of  fhe  curve. 

•  A  parlicularly  interesting  properly  of  fhe  proposed  approach  is  lhal  a  slochaslic  componenl 
may  easily  be  incorporated,  conlrary  to  fhe  continuous  curve  propagation  by  PDEs.  One  can 
indeed  add  a  random  perlurbalion  componenl  fo  fhe  deterministic  ODE  componenl,  such  as 
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a  standard  Brownian  motion  with  a  certain  variance,  in  order  to  make  the  evolution  more 
robust  to  noise  effects  which  usually  yield  a  convergence  to  a  local  extremum  of  the  objective 
functional. 

•  Finally,  we  incorporate  a  novel  global  regularizer  to  the  vertex  motion  equations  using  elec¬ 
trostatics  principles.  When  the  regularizer  is  the  sole  term  acting  on  the  polygon,  it  forces 
the  polygonal  edges  and  vertices  to  remain  apart  as  imposed  by  the  electric  field  induced 
at  a  vertex.  We  note  that  our  proposed  approach  is  to  be  distinguished  from  that  of  Bruck- 
stein  et.al  [20],  which,  using  a  discrete  form  of  curvature  to  account  for  only  local  geometry, 
evolves  a  polygon  in  the  absence  of  an  image  term  with  a  different  intended  application.  Our 
global,  rather  than  local,  regularizer  is  more  consistent  in  preserving  local  features  such  as 
corners. 


4.3  Active  Polygons 

The  dynamical  equation  of  an  active  contour,  typically  follows  the  construction  of  an  energy  func¬ 
tional  around  a  region  which  is  subsequently  minimized  by  a  gradient  descent  flow.  Our  goal  is  fo 
design  flows  to  move  a  polygon  by  its  relatively  small  number  of  vertices  rather  than  a  continuous 
active  contour. 

Let  us  consider  a  contour  energy  that  consists  of  an  integral  inside  and  outside  the  active  contour, 
where  the  integrand  f{x,y)  consists  of  the  function  /  :  — )■  M.  To  explicitly  invoke  a  contour 

C  :  [a,  6]  C  M  — )■  around  some  region  i?  C  we  use  the  divergence  theorem  to  write  an 
integral  over  the  interior  of  a  curve  as  a  contour  integral 

EiC)=  [[  f{x,y)dxdy=  i  {F,N)ds, 

JJR  JC=aR 

where  N  denotes  the  outward  unit  normal  to  C ,  ds  the  Euclidean  arclength  element,  and  where 
F  =  {F^,  is  chosen  so  that  V  •  F  =  /.  In  what  follows  we  will  let  p  G  [a,  b]  denote  a  fixed 
parameterization  of  the  curve  where  C  (a)  =  C  (b).  We  will  indicate  by  v  any  variable  whose 
variation  affects  the  geometry  of  the  curve.  In  the  case  of  a  polygon,  v  will  denote  a  Cartesian 
coordinate  of  any  vertex.  In  the  case  of  a  smooth  curve,  v  will  denote  a  curvilinear  coordinate 
which  varies  along  the  normal  direction  at  any  point  on  the  curve  but  remains  constant  along  the 
curve  itself.  The  key  point,  in  either  case,  is  that  v  and  p  constitute  independent  variables.  We  now 
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rewrite  the  eontour  integral  in  terms  of  the  parameter  p. 


E(C)=  /  {F,N)\\Cp\\dp=  /  {F,JCp)dp, 


where  J  = 


0  1 


and  fV||Cp||=JCp.  In  order  to  determine  the  gradient  flow  assoeiated 


L-1  0. 

with  E,  we  eompute  the  derivative  with  respeet  to  v.  After  some  manipulations,  Ey  can  be  obtained 
as  [150, 153] 

[■b 

(4.1) 


Ev{C)  =  f  f{Cy,JCp)dp. 
J  a 


In  the  case  of  a  smooth  curve  C ,  the  variable  v  denotes,  at  each  point  on  the  curve,  a  coordinate 
which  varies  in  the  normal  direction  to  the  curve.  Thus,  if  we  are  considering  geometric  evolutions 
of  the  curve,  we  see  from  the  final  expression  of  Ey,  that  the  gradient  flow  for  C  with  respect 
to  E  is  given  by 


dC 


fN. 


(4.2) 


Proceeding  along  a  similar  but  slightly  modified  line  of  thought,  we  consider  a  closed  polygon 
V  as  the  contour  C ,  and  with  a  fixed  number  of  vertices  {V  i, ..,  V  =  {{xi,yi),i  =  1, . . .  n}. 
We  may  parameterize  C  by  p  G  [0,  ra]  as 


C{p,V)=L{p-{p\,V^^,V  LpJ+i)  (4.3) 

where  {p\  denotes  the  largest  integer  which  is  not  greater  than  p,  and  where  L{t,  A  ,  B  )  =  (1  — 
t)A  +  tB  parameterizes  between  0  to  1  the  line  from  A  to  B  with  constant  speed,  (A  and  B 
denote  the  end  points  of  a  polygon  edge).  Note  that  the  indices  of  V  should  be  interpreted  as 
modulo  n  so  that  V  o  and  V  „  denote  the  same  vertex  (recall  C  is  a  closed  curve).  Finally,  note 
that  C  pis  defined  almost  everywhere  (where  p  /  LpJ  )  by 

C'p(p,V)  =  (4-4) 

Proposition  3  Using  the  above  parameterization  C  p{p,V  )  in  Eq.  (4.1),  we  obtain  the  first  varia¬ 
tion  of  the  energy  functional  E  as 

pn 

Ev=  f{Hp-\E\,V[j,\,V[p\+i)){Cv,J{Vipi+i-Vi^^))dp,  (4.5) 

Jo 

and  the  minimization  of  E  is  achieved  by  a  gradient  descent  flow  given  by  a  set  of  ODEs  for  each 
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vertex  V^,  k  =  1,  as  (see  Appendix B) 


=  Pf{Lip,  V  k-uVk))  dp  Nk,k_i  +  (1  -  P)f{L{p,  Vk,  V;fc+i))  dp  Nk+i,k,  (4.6) 

where  N^^k-i  (resp.  N^+i^k)  denotes  the  outward  unit  normal  of  edge  {V  k-i  —  Vk)  (resp.  (Vk- 

Vk+i)). 


Written  for  each  vertex,  these  equations  are  a  set  of  coupled  ordinary  differential  equations  to 
be  simultaneously  solved.  Intuitively,  any  of  these  equations  may  concisely  be  re-written  as 


dVk 

dt 


fk,k-ld^k,k-l  + 


(4.7) 


integrates  the  information,  which  is  obtained  from  image  values  along  two  adjacent  edges  {Vk-i  — 
Vk),  and  {Vk  —  Vk+i),  combined  with  the  global  image  statistics  (depending  on  the  function  /), 
into  two  overall  speed  functions  fk,k-i,  and  fk+i,k  iri  order  to  move  the  vertex  Vk  ■  This  system  of 
ODEs  hence  effectively  implements  a  coupled  motion  of  all  the  vertices  of  a  polygon.  In  addition, 
this  integration  procedure,  provides  improved  robustness  to  noise.  Another  advantage  of  the  flow 
in  Eq.(4.6)  is  the  reduction  of  the  dimension  of  contour  propagation  problem  from  a  theoretically 
infinite  one  to  roughly  on  the  order  of  3  —  30  vertex  points,  depending  on  the  complexity  of  an  object 
boundary.  This  also  clearly  highlights  the  differences  mentioned  earlier  between  active  polygons 
and  marker  particle  methods.  One  may  note  the  significantly  reduced  number  of  well  separated 
vertices  which  have  to  be  propagated.  In  addition,  the  motion  of  each  vertex  being  based  on  a 
weighted  combination  of  the  unit  normals  only  at  the  polygon’s  edge  points,  implies  that  there  is  no 
need  for  a  unit  normal  to  be  defined  at  a  vertex.  As  a  result,  the  polygonal  contour  formed  by  these 
vertex  locations  (Vj, ...,  Vjj)  need  not  be  differentiable.  A  generic  image-based  term,  indeed  a  speed 
function  /,  can  be  used  in  the  derived  ODEs  describing  the  motion  of  our  active  polygons,  hence 
making  them  quite  flexible. 

Eor  illustration,  we  apply  our  active  polygon  model  in  Eq.(4.6)  to  an  image  using  simple  mean 
statistics  for  the  foreground  and  the  background.  Eor  a  given  image,  with  a  simple  triangle  shape 
as  the  target,  the  propagation  of  an  active  polygon  with  different  initializations  toward  the  target 
shape  are  shown  in  each  row  of  Eig.  4.2.  The  edges  of  the  active  polygon  align  themselves  along 
the  edges  of  the  target  object,  with  the  influence  of  the  weights  determined  by  the  data  term,  i.e.  the 
two  weighted  integrals  along  two  neighboring  sides  at  each  vertex.  The  result  of  our  algorithm  for 
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Figure  4.2:  We  demonstrate  flow  (4.6)  for  two  different  initial  active  polygons. 


two  different  initializations,  given  at  the  right  of  each  row,  is  indeed  the  location  of  3  vertices  of  the 
target  triangle.  Implementation  issues  along  with  the  initialization  are  discussed  in  Section  4.6.1. 


4.4  An  Information-Theoretic  Criterion 

It  is  commonly  assumed  in  region-based  active  contours  that  an  image  is  piecewise  constant.  Sev¬ 
eral  techniques  based  on  utilizing  the  first  and  second  order  statistics  have  been  proposed.  While 
these  techniques  are  quite  adequate  for  Gaussian  data,  they  are  highly  insufficient  to  capture  the  un¬ 
derlying  information  in  non-Gaussian  data  which  include  almost  all  textures.  Towards  that  end,  we 
use  an  information  theoretic  measure  which  not  only  captures  the  higher  order  information  of  non- 
Gaussian  data,  but  provides  a  probabilistic  disparity  measure  among  N  data  populations.  Hence,  we 
consider  a  decision  problem  with  N  classes  ci, . . .  ,cn  with  prior  probabilities,  a  =  (ai, . . . ,  qn) 
such  that  —  1-  Based  on  the  Shannon  entropy  H,  the  so-called  generalized  Jensen-Shannon 

divergence  [90]  among  N  probability  densities  is 

(N  \  N 

(4.8) 

i=l  )  i=l 

where  Pi(^)  denotes  the  probability  density  of  the  class  in  a  region.  One  of  the  major  features  of 
the  Jensen-Shannon  divergence  is  that  different  weights  ai  may  be  assigned  to  the  relevant  distribu- 
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tions  according  to  their  importance. 

Our  ability  to  better  capture  the  underlying  statistics  of  regions  using  the  entropy-based  symmet¬ 
ric  divergence  measure,  suggests  that  this  may  play  a  key  role  in  constructing  an  energy  functional 
whose  optimization  would  yield  a  polygonal  flow  to  well  delineate  differently  textured  regions.  To 
proceed,  and  for  the  sake  of  efficiency,  we  resort  to  a  fast  numerical  estimation  of  the  densities  which 
in  turn  facilitates  the  estimation  of  the  entropies  and  hence  of  the  divergence  measure.  We  chose  to 
adopt  a  first  order  approximation  of  a  density  which  achieves  the  maximum  entropy  solution,  and 
which,  in  turn,  is  used  in  approximating  the  entropy  expression  as  proposed  by  Hyvarinen  in  [64]. 

On  a  region  delineated  by  an  active  contour,  we  assume  a  scalar  random  variable  (r.v.)  ^  on  a 
given  set  of  intensity  values  S,  ^  :  S  C  M  — )■  ,  and  the  available  information  on  the  density  of 
the  r.v.  ^  is  given  by 

J^PiOGjiOd^  =  Uj,  for  j  =  (4.9) 

where  the  estimates  Uj  are  the  expectations  E{Gj{^  )}  of  m  known  independent  functions 
of  ^  .  Note  that  there  is  no  model  assumption  for  the  random  variable  ^ ,  however,  the  distribution 
which  has  the  maximum  entropy  and  which  is  also  compatible  with  the  measurements  in  Eq.  (4.9) 
is  sought  [36,  69].  This  problem  has  been  widely  studied,  and  the  form  of  the  maximum  entropy 
distribution  has  been  derived  in  [36, 69].  The  density  PoiO  which  satisfies  the  constraints  (4.9)  and 
has  maximum  entropy  among  all  such  densities  is  of  the  formpo(0  =  where  B 

and  bj  are  constants  that  are  determined  from  the  constraints.  A  first  order  approximation  for  this 
maximum  entropy  density  denoted  by  p(^),  is  given  by  [64] 

(4.10) 

where  Uj  =  E{Gj{^)},  and  (/)(^)  is  the  standard  Gaussian  density  (/)(^)  =  exp(— ^^/2)/-v/^. 
(£'{•}  is  the  expectation  of  a  random  variable).  In  practice,  the  u/s  are  estimated  as  the  corre¬ 
sponding  sample  averages  of  the  Gj{^),  i.e.  we  compute  measurements  in  a  region  R  of  an  image 
function  7  :  — )■  M  by  Uj  =  '12{x  y)eRGjid{x,y)),  where  |7i!|  is  the  number  of  pixels  in 

R.  A  simple  approximation  of  the  entropy  functional  is  subsequently  found  upon  substituting  the 
approximate  density  p{^)  in 

r  1  ™ 

HiPiO)  =  -  JiO^ogpiOd^  ^H{v)  -  2  ^4.11) 


p{o  =  <p{o  1+E 


i=i 
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where  H^u)  =  ^(1  +  log(27r))  is  the  entropy  of  a  standardized  Gaussian  variable  (see  [64]  for 
details).  This  result  implies  that  a  first  order  approximate  maximum  entropy  of  a  given  distribution 
is  found  by  ealeulating  how  far  away  it  is  from  that  of  the  standard  Gaussian  density. 

4.4.1  Evolution  of  A  Single  Active  Contour 

With  a  single  eontour,  the  image  domain  C  is  split  into  two  regions,  namely  a  region  inside 
the  eontour,  eall  it  R^,  and  a  region  outside  the  eontour,  0\Ru-  Exploiting  the  approximation  in 
Eq.  (4.11),  we  define  an  energy  funefional  whose  opfimizafion  yields  fhe  evolution  of  our  aefive 
eontour  based  on  a  divergenee  measure  defined  in  Eq.  (4.8).  The  new  energy  funefional  for  fwo 
regions,  is 


JSa,2  =  iT(aiPl(0  +«2P2(0)  -  ^(Pl(0)  -  «2 

=  -  JjaiPliO  +«2P2(^))log(aiPl(0  +  «2P2(0)c^^ 

+  (31  ^og(pi(^))d^  +  a2  J^P2{0  log(p2{0)d^, 


(4.12) 


Here  fhe  expression  H{ai  pi{^)  +  (32  P2{0)  similarly  to  H{pi{^))  be  approximated  by  sub- 
sfiluling  fhe  approximate  density  expression  in  Eq.(4.10)  fo  yield 


Hiaipii^,)  +a2P2iO)  ~  - 


i=i 


A  III  III 

«i  +  a2  (/((^[l  +  Y^VjGj{0] 

i=i  i=i 

(m  m  \ 

«1  +  YlujGj{0]  +  a2  +  1 

-  Uj  +  (32  Vj)  Gj(oj 

'  log  <PiO  +^(ai  Uj  +  (32  Vj) 


^  m 

H{i2)  -  + 


(32  Vj  ) 


(4.13) 


i=i 


This  shows  fhaf  fhe  same  firsl  order  enfropy  approximation  holds  for  fhe  sum  of  densities.  Denof- 
ing  measuremenfs  inside  fhe  eontour  as  Uj,  and  fhose  oufside  as  Vj,  for  j  =  1,  ..,m,  and  using 
fhe  approximate  enfropy  expressions  (4.11),  (4.13)  in  fhe  energy  funefional  Eq.(4.12),  if  may  be 
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approximated  as 


1  m  / 

H{v)  -  -  ^(ai  Uj  +  a2  vjf  -  ai  H{v) 

i=i  V 

^  ra 

-  ^  (-(ai  Uj  +  (32  +  (3l  +  (32  Uj) 

i=i 

^  m 

-  ai  (32  XlK'  “ 

i=i 


(32 


1 

2 


E 


(4J4) 


by  noting  that  ai  +  (32  =  1- 

Note  that  on  aeeount  of  the  higher  order  nature  of  the  eoeffieients  Uj  and  Vj  (i.e.  higher  order 
than  first  and  seeond  moments),  the  proposed  energy  funetional  subsumes  the  previously  proposed 
teehniques  based  on  the  first  and  seeond  order  statisties.  Those  may  in  faet  be  shown  to  be  partieular 
eases  of  the  above.  This  may  also  be  justified  by  the  faet  that  non-Gaussian  densities  may  be 
expanded  and  that  their  higher  order  eumulants  refieet  the  degree  (or  the  laek  thereof)  of  skewness  or 
kurtosis  of  a  density  relative  to  a  normal  [71],  henee  G(^)  =  and  G(^)  =  would  be  partieular 
ehoiees  for  the  measurement  funetions  Gj.  Other  ehoiees  of  these  funetions;  Gi  (^)  =  as  an 

odd  funetion  to  measure  asymmetry  (analogous  to  the  third  moment  as  a  measure  of  skewness),  and 
G2{0  —  1^1 5  oi"  as  ehoiees  of  even  funetions  (analogous  to  the  fourth  moment  as  a  measure 

of  sparsity,  bimodality,  relative  to  a  Gaussian  density),  are  given  in  [64,  65].  The  idea  behind  using 
expeetations  of  odd  and  even  funetions  of  the  data  is  thus  similar  to  attempts  to  eharaeterize  a  density 
by  the  first  few  moments  (usually  order  is  less  than  4).  The  ehoiee  of  priors  ai  for  eaeh  density  is 
explained  next: 


-  If  we  assign  eonstant  priors  ai  and  a2  (e.g.  equal  priors  ai  =  a2  =  0.5)  for  both  of  the 
regions,  the  eonstant  (3i(32/2  multiplying  the  energy  funetional  in  Eq.  (4.14)  has  no  effeet  in 
terms  of  its  first  variation,  and  the  gradient  deseent  flow  of  the  aetive  eontour  that  minimizes 
the  energy  funetional  E  =  —  JSa,2  =  —^u,i  (32  ^]Liiuj  —  Vj)‘^  may  be  obtained  by  taking 
its  first  variation  as 


—  =  VJSa,2  =  ^{Uj  -  Vj)(Vuj  -  Vvj),  (4.15) 

i=i 

where  a  eontinuous  aetive  eontour  C  is  utilized.  As  we  have  noted  in  the  previous  seetion, 
the  resulting  flow  ean  be  direetly  applied  to  an  aetive  polygon  instead  of  a  eontinuous  aetive 
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contour.  A  measurement  or  a  constraint  on  a  region  distribution;  Uj  (resp.  Vj)  for  region  Ru 
(resp.  Ry),  is  given  by 


Gj(I(x,y))dxdy  _  Gj{I{x,y))dxdy 


\Rif 


\Ri, 


(4.16) 


with  \Ru\  =  dxdy,  dxdy,  and  for  j  =  l,...m  different  constraints.  Then 

the  partial  variation  of  Uj  and  Vj  in  Eq.(4.16)  w.r.t  C  is  given  by 


VcUj  = 


Gj (I)  Uj 

\Ru\ 


5  ^  c'^j  - 


\R.\ 


Nu, 


(4.17) 


where  Ny  denotes  the  outward  unit  normal  of  C  (region  Ry).  Note  that  the  outward  unit 
normal  for  the  boundary  of  the  outer  region  is  —Ny. 

The  contour  evolution  is  found  by  substituting  Eqs.(4.17)  into  the  gradient  descent  flow 
Eq.(4.15): 


dCu 

dt 


m 


fNu,  where  /  =  >  (r/j  -  vj) 


i=i 


Gj{I{x,y))  -  Uj  Gj(I(x,y))  - 


\Ry 


+ 


\Rv 


(4.18) 


We  note  that  this  is  a  generalized  form  of  the  data  term  of  the  flow  proposed  by  Yezzi,  Tsai, 
and  Willsky  [149, 150].  In  the  above  equation,  the  speed  /  of  the  contour  along  its  normal 
direction,  can  be  directly  used  in  the  active  polygon  evolution  equation  derived  in  Eq.  (4.6), 
and  restated  in  Eq.  (4.7). 


-  On  the  other  hand,  in  Eq.  (4.14)  we  may  assign  variable  weights  ai,  and  a2  that  depend  on 
the  regions.  An  intuitive  choice  would  be  to  pick  the  ratio  of  the  area  of  each  region  to  the 
total  area  of  the  image  domain,  say  A  =  \Ry  \  +  \Rv\.  This  choice  in  practice  implies  taking 
into  account  the  number  of  pixels  in  each  region,  and  would  thus  lead  to  a  measure  that  is 
normalized  with  respect  to  the  areas  of  the  regions.  Eetting  ai  =  a2  =  '"a’ 
have  to  take  into  account  in  the  derivation  of  the  gradient  descent  flow,  the  first  variations  of 
the  coefficient  terms  as  well: 


dC 

~dt 


VJS^ 


a, 2 


E 

i=i 

1  \Rv 


1  \Ru\\R„ 


2  A2 


{2{uj  —  Vj){Vuj  —  Vvj))  Nu 

1  \R. 


2  A2 


(Ui  -  vA  Nu  -  - 


2  A2 


(ui  -  vA  Nu- 
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+ 


(4.19) 


Figure  4.3:  A  continuous  contour  with  a  small  regularization  (top  row),  fails  to  capture  as  a  whole 
an  object  with  a  synthetic  texture  of  vertical  stripes  ,  whereas  with  a  large  regularization  (middle 
row),  rounds  off  the  boundaries,  and  continues  to  shrink  without  sticking  to  the  data.  A  polygonal 
contour  (bottom  row)  correctly  segments  the  textured  region. 


After  some  manipulations  (given  in  Appendix  C),  this  gradient  descent  flow  can  be  written  as 
dC  1 

—  =  fNu,  where  /  =  ^  '^{uj  -  Vj)  -  uj)  +  (Gj(I)  -  vj))  (4.20) 

j=i 

whose  energy  functional  is  indeed  a  generalized  form  of  the  external  energy  functional  that 
Chan  and  Vese  proposed  in  [29]: 

Ecv  =  ^f2  [  (Gjil)  -  ujfdxdy  +  f  (Gjil)  -  Vjfdxdy.  (4.21) 

We  also  note  here  that  in  the  presence  of  severe  noise,  as  illustrated  in  Fig.  4. 1  for  a  uniform 
noise,  the  continuous  contour  ends  up  encircling  small  noisy  regions,  and  its  length  grows.  To 
overcome  such  effects,  a  penalty  on  the  length  of  an  active  contour  is  added  to  its  energy  func¬ 
tional  E  =  JJj^  fdxdy  +  a  ds,  where  o;  >  0,  and  s  is  the  arclength  of  C  (also  explained  in 
Section  2.3.3).  This  second  term  brings  a  problematic  trade-off  which  we  avoid  in  active  polygon 
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Figure  4.4:  An  active  contour  (first  row)  fails  to  capture  synthetic  texture  of  vertical  stripes  even 
after  a  Gabor  filtering,  whereas  the  active  polygon  (second  row)  captures  the  outline  of  the  textured 
region. 

framework  as  will  be  explained  next. 

The  advantage  of  propagating  a  polygon  instead  of  a  continuous  curve,  is  demonstrated  in 
Fig.4.3,  where  a  synthetic  texture  with  vertical  stripes  is  given.  Obviously,  thresholding  for  seg¬ 
mentation  will  not  work  since  one  of  the  stripes  has  the  same  color  as  the  background.  A  contin¬ 
uous  active  contour  propagation  (first  row  of  Fig.4.3)  fails  with  a  small  regularization,  i.e.  small  a 
mentioned  in  the  previous  paragraph,  because  the  curve  dips  down  into  the  gaps  of  the  stripes  with 
the  same  color  as  the  background,  and  the  individual  bars  are  captured.  With  an  increased  amount 
of  regularization  (higher  a)  on  the  continuous  curve  (second  row  of  Fig.4.3),  the  curve  remains  in¬ 
tact,  however  comers  are  rounded,  and  if  one  maintains  the  evolution,  the  important  features  will  be 
“missed”.  We  have  not  been  successful  at  obtaining  the  appropriate  tradeoff  between  the  data  term 
and  the  regularizer  in  the  continuous  case  to  yield  a  satisfactory  result.  Active  polygon  propagation 
in  Eq.(4.6)  with  /  in  (4.20)  (third  row  of  Fig.4.3),  accurately  and  consistently  captures  the  target 
shape.  We  also  show  in  Fig.4.4,  the  Gabor-filtered  version  of  the  same  synthetic  textured  object, 
(tuned  for  vertical  orientation).  It,  however,  can  be  observed  that  the  active  contour  still  fails  to 
operate  on  the  filtered  texture,  whereas  the  active  polygon  again  successfully  captures  the  outline  of 
the  target  region. 
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4.4.2  Evolution  of  Multiple  Active  Contours 


A  nice  property  of  the  Jensen-Shannon  divergence  measure  is  that  it  may  be  generalized  as  a  prob¬ 
abilistic  difference  measure  among  any  finite  number  of  probability  densities.  Coupled  propagation 
equations  for  multiple  active  contours  which  delineate  different  regions  on  an  image  domain  can 
also  be  obtained  as  the  gradient  ascent  flow  of  the  JS  divergence  measure  among  the  densities  of 
those  regions  (note  that  we  are  trying  to  maximize  the  distance  among  the  densities  of  regions). 
There  is  flexibility  in  the  placement  of  the  contours,  which  may  delineate  distinct  or  overlapping 
regions  without  difficulty.  We  again  derive  the  gradient  flows  generally  with  respect  to  continuous 
contours.  Figure  4.5  depicts  two  active  contours  whose  inner  regions  are  denoted  by  R^,  and  Ry, 


Figure  4.5:  Ternary  image  regions. 


and  their  common  exterior  by  Ry,.  The  measurements,  i.e.  statistics,  in  these  respective  regions 
are  denoted  hy  j  =  1, ...  ,m  for  m  different  measurements,  with  respective  prior  prob¬ 
abilities  di,  (32,  «3  :  ai  +  a2  +  =  1.  The  energy  functional  for  three  densities  can  be  written 

as 


JS 


a, 3 


i=l 


^  m 

H{v)  -  -  ^(ai  Uj  +  a2  vj  + 


as  Wj) 


i=i 

-ai  \  H(u)  - 


I  -“2 


-as  H{v) 


^  m 

-  ^  {ai  a2{uj  -  Vjf  +  ai  as{uj  -  Wjf  +  a2  as{vj  -  Wjf) 


i=i 


(4.22) 


This  form  of  energy  functional  may  easily  be  extrapolated  to  three  active  contours,  thus  four 
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regions  with  statistics  Uj,Vj,Wj,  Zj, 


m 


E 


0.1  02{uj  —  Vj)"^  +  ai  az{uj  —  Wj)"^  +  ai  ai{uj  —  ZjY 


+  (32  03  {Vj  -  WjY  +a2  (34  (wj  -  ZjY  +  (34  (Wj  -ZjY-  (4.23) 


The  prior  probabilities  of  each  density  ai, . . . ,  (34  may  be  selected  in  a  variety  of  ways  with  the 
simplest  selection  being  the  equal  priors,  i.e.  =  0.25,  *  =  1, . . . ,  4.  It  is  now  straightforward  to 
approximate  the  JS  functional  J S a, n  for  N  regions  on  the  image  domain.  It  may  be  observed  that 
the  first-order  approximations  to  both  the  densities  and  the  corresponding  entropies  of  the  regions, 
lead  to  an  overall  measure  that  computes  a  weighted  sum  of  the  divergence  measures  (i.e.  distances 
between  their  statistics)  between  all  pairwise  combinations  of  the  regions. 

Active  contour  evolutions  for  three  regions  using  means  and  variances  as  statistics  of  each  region 
were  derived  from  a  totally  different  perspective  by  Yezzi,  Tsai,  and  Willsky  [149].  Their  energy 
functional  was  based  on  a  geometric  notion  which  maximized  the  area  of  a  triangle  formed  by  the 
statistics  of  the  three  regions.  Our  energy  functional  on  the  other  hand,  is  information-theoretic  in 
nature,  and  evaluates  the  distance  among  probability  densities. 

The  ternary  case  entails  the  derivation  of  a  gradient  flow  for  each  of  the  two  active  contours 
C  u,  and  C  y.  Taking  the  first  variation  of  the  energy  functional  in  Eq.(4.22)  w.r.t.  the  contour  C„, 
yields 

m 

Vc„^*S'a,3  =  '^[oi  02iUj-Vj)+ai  03{Uj-Wj)]Vcu'Oj-[oi  a3(Uj-Wj)  +  a2  OsiVj-WjYVcuWj 
i=i 

(4.24) 

where  we  used  the  fact  that  =  0,  Vj,  since  Vj  are  the  statistics  inside  the  contour  C y  which 

do  not  depend  on  the  contour  C y. 

The  partial  variation  of  Wj ’s  requires  more  attention  than  that  of  Uj ’s  and  Vj ’s  since  the  statistic 
Wj  is  calculated  over  the  common  exterior  of  both  contours  whose  boundary  may  not  be  smooth 
when  the  two  contours  overlap.  We  exploit  a  a  similar  strategy  given  in  [149]  to  express  a  statistic 
in  this  third  region  using  characteristic  functions  Xu^Xv  over  the  regions  Ry,  Ry  as 

In\R^  Gj{I){l  -  Xv)dxdy 

qjj  ■  ^  _ i _ 

^  -  Xv)dxdy 

where  the  denominator  may  be  renamed  as  \Rw\-  The  variation  of  Wj  w.r.t.  C  y  can  hence  be 
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obtained  as 


=  - 


\Rw\ 


(1  -  Xv)Nu- 


(4.25) 


The  gradient  descent  flow  of  the  contour  which  is  equal  to  Vcu  JSa,z  is  thus  obtained  as 
'^here 

V^r  /  \  !  \tG  j{I{x,y))  —  Uj 

“2(%'  -  Vj)  +  ai  as{uj  -  Wj)\^ - — — - 2. 

i=i  ' 

+  [ai  asiuj  -  Wj)  +  a2  az{vj  -  — ^(1  -  Xv)-  (4.26) 

\-^W  \ 

By  similar  arguments,  one  can  write  Wi  =  where  the  denominator  is 

^  fn\Ry(^-Xu)dxdy 

equal  to  \R-w\,  the  variation  of  Wj  w.r.t  C y  can  be  obtained  as 


-  Xu)Ny. 


(4.27) 


Then  taking  the  first  variation  of  the  energy  functional  in  Eq.(4.22)  w.r.t.  the  contour  C  y,  and 
noting  that  Vc„Uj  =  0,  Vj,  leads  to  the  gradient  flow  of  the  contour  C  y,  ^  =  fyNy,  where 


fv  =  '^[-ai  a2iuj  -  Vj)  +  a2  azivj  -  Wj] 

i=i 


Gj{I{x,y))  -  Vj 

'  \R„\ 


r  /  N  /  \iGi{I{x,y))  —  Wj  .  . 

+  [ai  ari{uj  -  Wj)  +  a2  aji{vj  -  Wj)\^ - — — - ^(1  -  Xu)- 

■tlin 


(4.28) 


The  two  speed  functions  /„  and  fy  when  inserted  into  Eq.(4.6)  for  two  separate  polygon  propaga¬ 
tions,  C  u  and  C y,  result  in  ternary  polygonal  flows. 

We  illustrate  in  Eigure  4.6,  the  ternary  case  for  polygon  propagations  with  /„  and  fy  given  in 
Eq.(4.26)  and  Eq.(4.28).  Although,  the  polygonal  contours  and  C  y  evolve  separately,  their 
motion  is  coupled  through  the  variables  in  the  evolution  equations  that  depend  on  each  of  the  three 
regions.  Two  contours  move  in  such  a  way  to  maximize  the  approximate  Jensen-Shannon  diver¬ 
gence  among  densities  of  the  three  regions,  namely  Ry.  inside  the  contour  Ry:  inside  the 
contour  C  y ;  Ry, :  the  complement  of  Ry\^Ry.  Only  means  are  used  as  the  separating  statistics, 
(i.e.,  j  =  1,  G{x)  =  x),  and  the  resulting  two  polygons  are  shown  in  Eigure  4.6  (right).  Hence,  the 
gain  is  again  two-fold:  segmentation  of  the  targets,  and  their  description  in  terms  of  a  handful  of 
vertices  are  both  achieved. 
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Figure  4.6:  Ternary  flows  using  image  forces  with  G{x)  =  x,  are  used  to  segment  this  simple 
ternary  image  corrupted  by  Gaussian  noise.  Output  of  the  algorithm  shown  on  the  right  provides 
only  5  vertex  locations  for  the  first  object,  and  3  vertex  locations  for  the  second  object. 


4.4.3  Active  Contours  for  Vector- Valued  Images 

Up  to  now,  we  have  derived  the  evolutions  for  grayscale  (scalar  intensity)  images.  These  results 
may  be  extended  to  vector-valued  imagery  (e.g.  RGB  (Red,Green,Blue),  or  other  multi-spectral) 
with  rich  number  of  different  channels.  The  vector-valued  image  I  ^  ,  can  be  expressed 

in  terms  of  its  component  functions  :  I  =  {I{x,  y,l),  I{x,y,2), . . .  ,  I{x,  y,  rich)-  The  two  speed 
functions  for  the  ternary  case,  may  for  instance,  be  computed  as 


rich  m 


fu  =  a2iuj{l)  -  Vj{l))  +  ai  a^iujU)  -  Wj{l))] 

1=1  j=i 


Gj{I{x,y,l))  -  Ujil) 

\Rn\ 


+[ai  az{uj{l)  -  Wj{l))  +  a2  az{vj{l)  -  -»;j(0)]  ^ , — ^^^(1  -  Xv)}, 


fv  =  J2^{[-ai  a2{uj{l)  -  Vj{l))  +  a2  aaivjil)  -  Wj{l))]^ - — - 

1=1  j=l  ' 

-f [ai  aziujil)  -  Wj{l))  +  a2  (4.29) 

-^7/1 


4.5  A  Global  Polygon  Regularizer 

The  flow  of  an  active  polygon  may,  under  the  sole  influence  of  a  data  term,  become  undefined 
(degenerate)  under  a  variety  of  scenarios,  e.g.  when  a  vertex  becomes  infinitesimally  close  to  a 
non-neighbor  edge  of  the  polygon,  or  when  two  vertices  or  two  edges  come  infinitesimally  close  to 
each  other  at  some  point.  We  show  such  an  example  in  Fig.4.7,  in  which  propagation  of  an  initial 
polygon  with  5  vertices  becomes  degenerate,  where  the  vertex  V 4  approaches  (then  crosses)  the 
edge  between  vertices  V  2  and  V3 .  Asa  solution  to  overcome  this  problem,  we  introduce  a  natural 
regularizing  term  well  adapted  to  an  evolving  polygon.  Note  that  the  regularizer  we  are  about  to 
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Figure  4.7:  We  demonstrate  that  flow  (4.6)  may  become  degenerate  without  an  additional  constraint 
on  the  motion  of  vertices. 

propose  is  highly  adapted  to  propagating  polygons,  and  would  not  be  computationally  feasible  for 
smooth  continuous  contours.  This  regularization  is  accomplished  by  first  viewing  each  edge  of  the 
polygon  as  a  finite  line  charge  of  uniform  charge  density.  An  electric  field  E ,  generated  by  a  finite 
uniform  line  charge,  exerted  on  a  point  charge,  is  directed  away  from  the  test  charge  (we  assume  all 
polygon  edges  have  same  charge  polarity). 

4.5.1  Electric  force  by  a  line  charge: 

Given  a  line  charge  or  a  rod  of  positive  charge  that  extends  from  a  generic  point  a  G  to  6  G  , 
our  goal  is  to  calculate  the  electric  field  at  a  point  x'  G  .  Points  are  chosen  in  ,  and  using 
electrostatics  principles,  we  consequently  derive  these  fields  for  them.  The  line  charge  is  assumed 
to  be  made  up  of  differential  point  charges  dq.  We  need  to  compute  the  differential  electric  field 
dE  {x ')  exerted  at  a: '  by  a  charge  dq  at  location  x  =  a  +  t{b  —  a  )  which  is  on  the  rod  (depicted 
in  Figure  4.8).  As  given  by  Coulomb’s  law  [138],  dE  {x')  is  inversely  proportional  to  the  square 
of  the  Euclidean  distance  ||a:'  —  a:  |p  between  x  and  x',  and  its  direction  is  given  by  the  vector 
{x'  —  x)l\\x'  —  x\\,  i.e.,  dE  [x')  =  {x'  —  x)  l\\x'  —  x\\^dq.  We  assume  a  uniform  charge 
density  A  (a  constant  parameter)  along  the  rod,  hence  dq  =  \dx.  With  the  change  of  variable 
X  =  a  +t{b  —  a),  the  differential  amount  of  increase  in  dx  =  Ldt,  where  L  =  ||6  —  a  ||  is  the 
length  of  the  rod  extending  from  a  to  b. 

By  the  principle  of  superposition,  the  total  electric  field  E  acting  on  point  charge  at  x '  due 
to  line  charge  (a  ,  6 ),  can  be  obtained  by  integrating  the  fields  contributed  by  all  differential  point 
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dE(x’) 


Figure  4.8:  Calculation  of  electric  force  exerted  by  a  line  charge  (a  ,  6 )  at  a  point  x '  on  the  polygon. 


charges  on  the  rod  making  up  the  linear  charge  distribution, 

fb 
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=  kr  L 
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rdt. 


(4.30) 


Note  that  the  constant  A  is  combined  with  the  Coulomb’s  constant  into  a  constant  parameter  called 
kc-  After  some  manipulations  (see  Appendix  D),  the  total  electric  field  E  ab{x ')  can  be  written  in 
terms  of  the  two-vectors  x  a  =  x'  —  a  and  x  ^  =  x'  —  b  (depicted  in  Figure  4.8), 
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(4.31) 


As  shown  in  the  previous  section,  a  data  term  corresponding  to  a  given  vertex  consists  of  an 
integrated  force  along  its  two  neighboring  edges.  In  a  similar  way,  an  electric  field  exerted  af  a 
vertex  is  also  infegrafed  along  fwo  neighboring  edges.  Evolution  of  a  polygon  vertex  V  k,  due  fo 
fhe  fofal  elecfric  field  E  ^  integrated  along  fwo  neighbor  edges  {V  k)  and  {V  k,V  may 

fhen  be  written  as 


dVk 

dt 


Ek  = 


n  „1 

V  /  p  Ey  ,  y  XpV  k  +  0-  -  P)V  k-l)  dp 

Jo  3  '  3 

j  =  0 


j  X^k,ji^k  -  I 

n  „1 

+  XI  /  (1  +P^A:-rl)4'-  (4-32) 

J  0 

j  =  0 


j  ^  k,j  ^  k  +  1 


This  global,  rafher  fhan  local,  geomefric  dependence  makes  fhis  regularizer  very  differenf  from 


75 


Figure  4.9:  Electrostatic  regularizer  in  Eq.(4.32)  computes  the  electric  force  at  each  vertex. 

the  ones  used  in  the  literature.  The  use  of  a  global  regularizer  is  more  consistent  with  the  desire 
to  capture  local  features  that  would  otherwise  not  be  captured.  Thus,  the  additional  and  novel 
geometric  component  of  the  polygon  evolver,  which  induces  global  geometric  dependence,  provides 
a  regularization  which  both  avoids  the  flow  degeneracy  as  well  as  captures  sharp  comers  of  the 
target  shape  without  any  shrinking  or  smoothing  effects. 

We  demonstrate  the  polygon  regularization  capability  of  the  electric  force  flow  in  Eq.(4.32)  in 
Eigure  4.9  for  two  different  polygons.  The  polygonal  evolution  with  an  electrostatic  regularizer 
force  (plotted  at  each  vertex)  is  shown  at  several  snapshots  in  these  figures.  Note  how  the  force 
at  each  vertex  is  large  in  magnitude  initially,  and  push  a  vertex  away  from  the  other  edges  of  the 
polygon.  The  magnitude  of  the  electric  force  at  each  vertex  decreases  and  edges  remain  well  apart 
which  is  exactly  the  effect  we  expect  from  the  regularizer  force.  This  force  should  be  insignificant 
when  a  vertex  and  its  adjacent  edges  are  not  very  close  to  most  of  the  other  edges,  and  should 
become  influential,  even  dominate  when  the  vertex  or  its  adjacent  edges  are  very  close  to  other 
edges. 

The  addition  of  the  regularizer  term  (4.32)  to  the  motion  equation  of  a  vertex  obtained  in  Eq.(4.6) 
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Figure  4.10:  The  electric  forces  for  an  evolution  are  shown.  Note  the  forces  only  become  significant 
when  a  vertex  approaches  another  vertex  or  edge. 


leads  to  the  following  modified  verfex  flow 


Here,  a,  a  consfanf  paramefer  fo  weighl  fhe  influence  of  dafa  term,  and  fhe  elecfric  field  ferm,  is 
chosen  as  0.95  fhroughouf  fhe  evolutions.  The  reason  fhaf  we  pul  such  a  heavy  weighl  on  fhe  dafa 
ferm  is  fhaf  fhe  regularizer  only  kicks  in  very  powerfully  when  degeneracy  occurs,  and  if  lels  fhe 
dafa  term  govern  fhe  evolufion  of  fhe  polygon  during  mosl  of  fhe  evolufion  lime. 

In  Figure  4.10  we  demonslrale  fhe  use  of  flow  (4.33)  for  fhe  same  inilial  aclive  polygon  wilh 
5  vertices  which  has  been  shown  fo  resull  in  an  ill-posed  flow  in  fhe  previous  section  in  Figure 
4.7.  Here,  snapshols  of  only  fhe  active  polygon  on  fhe  friangle  shape  image  are  given  fo  belter 
appreciate  fhe  influence  of  fhe  regularizer,  and  we  show  fhe  elecfric  force  al  each  vertex  during  Ihis 
evolufion.  Note  fhaf  fhe  electric  force  at  a  vertex  becomes  significantly  large  when  the  vertex  is 
infinitesimally  close  to  another  edge  in  this  figure.  This  evenl  exaclly  keeps  fhe  polygon  simple 
during  fhe  evolution,  and  by  fhe  effecl  of  fhe  dafa  ferm,  fhe  polygon  converges  fo  fhe  largel  shape. 
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4.6  Results  and  Conclusions 


In  this  chapter,  we  have  presented  a  new  approach  for  image  segmentation  through  polygon  prop¬ 
agating  equations.  In  this  section,  by  way  of  numerical  experiments,  we  validate  the  effectiveness 
and  the  usefulness  of  our  technique. 

4.6.1  Implementation  Issues 

The  performance  of  most  active  contour  algorithms  depend  on  the  initial  conditions  [153].  A  partic¬ 
ularly  important  question  in  carrying  out  such  flows,  is  that  of  initializing  the  contour.  Specifically 
in  our  case,  we  need  to  specify  the  number  of  vertices,  and  their  placement  to  start  off  the  evolu¬ 
tion  of  the  polygon.  For  the  active  polygons,  though,  it  is  possible  to  circumvent  this  problem,  and 
speed  up  the  convergence,  by  a  simple  approach  which  helps  an  initial  active  polygon  adaptively 
adjust  to  the  number  of  vertices  required  for  the  description  of  the  target  shape.  Towards  that  end, 
we  may  initialize  a  very  coarse  polygon  with  a  small  number  of  vertices,  e.g.  a  big  rectangular  or 
circular  polygon  close  to  the  image  boundaries.  While  the  initial  polygon,  usually  with  a  very  few 
number  of  vertices,  is  propagating  with  both  the  image  force  and  the  regularizer  force,  new  vertices 
are  periodically  added  and  removed  affording  it  a  flexible  motion  towards  the  target  region.  A  nat¬ 
ural  criterion  to  remove  a  vertex  may  for  instance  be  based  on  the  angle  between  its  two  adjacent 
edges  being  close  to  either  0,  or  tt.  This  may  be  effected  for  a  given  vertex  V  k,  with  its  adjacent 
edges  A  =  {V  k-i  —  V  k),  and  B  =  {V  h+i  —  V  k)^^y  computing  the  following  inner  product 
AB  =\\A\\\\B  1 1  cos (Oab),  the  9ab  being  the  angle  between  the  two  vectors.  One  may  therefore 
check  if 

I  cos(0nB)|  =  Ill'll  I  «  1,  (4.34) 

to  determine  the  redundancy  of  a  vertex.  In  either  case,  the  vertex  V  k  may  simply  be  removed  from 
the  vertex  list,  to  finally  yield  a  polygon  properly  enclosing  the  shape.  During  a  vertex  addition 
period  of  an  evolution,  the  magnitude  of  our  image  force  along  each  edge  of  a  polygon,  i.e.  D  = 
Q\f{L{p,Vk,Vk+i  ))\dp,  is  computed,  and  a  new  vertex  is  added  to  a  middle  point  of  the  edge 
with  the  maximum  value  of  D.  The  intuition  here  is  that  the  edges  with  higher  image  speeds  are 
closer  to  image  structures  that  may  require  finer  details. 

For  implementational  purposes,  the  polygon  structure  P  is  a  two-vector  (2-D  vector),  and  the 
normals  needed  in  the  computation  of  Eq.(4.33)  are  computed  as  N^^k-i  =  {—P{k,  1)  +  P{k  — 
1,1),  P{k,0)  —  P{k  —  1,0))^.  We  also  maintain  a  two-dimensional  function  which  acts  as 
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an  indicator  function  to  show  which  pixel  is  inside  or  outside  a  polygon  during  the  evolution  to 
facilitate  computation  of  the  statistics  in  the  regions.  Statistical  calculations  may  be  carried  out 
fast  by  computing  the  change  in  the  position  of  each  polygon  edge,  and  appropriately  adding  or 
subtracting  the  statistics  computed  in  that  difference  region,  which  results  in  the  motion  of  the  edge, 
which  is  a  straight  line,  from  one  iteration  to  the  next.  Time  step  in  the  discretization  of  the  ODE 
is  chosen  large,  e.g.  5t  =  10,  which  indeed  helps  the  model  to  escape  various  local  minima  of  the 
objective  function.  For  a  typical  image  of  256  x  256  size,  the  relative  speed  of  the  active  polygon 
model  is  approximately  twice  that  of  an  associated  active  contour  model  using  the  same  data  term. 

As  we  mentioned,  we  may  also  add  a  random  perturbation  to  our  ODE  model,  a  zero-mean 
Gaussian  r.v.  with  a  variance  that  has  a  very  minor  effect  when  compared  to  the  change  in  the  image 
term.  This  trivial  perturbation  changes  the  path  of  each  evolution  very  slightly  although  the  results 
are  very  consistent  as  will  be  shown  next. 

4.6.2  Experimental  Results 

In  this  section,  we  demonstrate  texture  segmentation  examples  on  natural  texture  images.  Our 
active  polygon  propagation  model  is  in  Eq.(4.33)  with  speed  function  in  (4.20),  and  measurement 
functions  Gi(^)  =  G2(^)  =  G3  =  |^|,  G4  =  log(cos/i(rE)). 

In  the  first  example,  a  zebra  on  a  grassy  background  constitutes  of  mainly  two  textured  regions 
(Fig.4.1 1).  A  generic  rectangle,  i.e.  just  four  vertices,  is  initialized  on  the  zebra  image  which  is  quite 
challenging  in  terms  of  unsupervised  texture  segmentation.  Snapshots  from  the  polygon  propagation 
with  the  resulting  segmentation  in  Fig.4.11,  show  that  a  zebra  figure  is  very  nicely  capfured. 

Ofher  nafural  fexfure  examples  include  a  monarch  larvae  and  a  monarch  bullerfly  wifh  generic 
polygon  inifializafions,  a  circle  and  a  recfangle,  are  shown  in  Fig.4.12.  In  fhe  same  figure,  an- 
ofher  arbifrary  inifializafion  on  fhe  same  monarch  picfure,  show  fhaf  fhe  fargef  fexfured  body  of  fhe 
monarch  is  capfured  in  bofh  cases.  In  Fig.4.13,  a  fish  whose  body  has  a  fexfure  of  sfripes  is  cap¬ 
fured  by  an  active  polygon.  Similarly,  a  sea  sfar  on  a  fexfured  rocky  ferrain,  (Fig.4.14),  a  cheefah 
figure  (Fig.4.15  leff),  anofher  cheefah  in  bushes  (Fig.4.15  righf),  and  a  chunk  of  crysfal  (Fig.4.16) 
are  shown  fo  demonsfrafe  fexfure  capfuring  capabilifies  of  fhe  acfive  polygon  model  fogefher  wifh 
fhe  Jensen-Shannon  criferion. 

One  of  our  goals  was  fo  apply  our  active  polygon  model  fo  capfure  man-made  objecf  shapes.  A 
group  of  real  airplane  image  experimenfs  on  each  of  which  a  generic  circular  polygon  was  inifial- 
ized,  are  shown  in  Fig.4.17.  Our  model  successfully  capfures  fhe  plane  shape  in  terms  of  a  polygon. 
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Figure  4.11:  A  zebra  figure  is  captured  by  the  active  polygon  model.  A  generic  rectangular  active 
polygon  close  to  image  boundaries  is  initialized. 


80 


Figure  4.12:  A  monarch  larvae  on  a  leaf  is  captured  by  an  active  polygon  (left).  Monarch  butterfly 
is  captured  in  the  two  other  columns  with  very  different  initalizations. 
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Figure  4.13:  A  fish  with  a  striped  texture  is  captured. 
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Figure  4.15:  Left:  Cheetah  figure  is  captured  by  the  active  polygon.  Right:  A  cheetah  in  bushes  is 
captured  with  an  active  polygon. 
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which  may  be  essential  as  an  input  to  object  recognition  algorithms.  An  image  of  another  man-made 
object,  a  submarine,  is  shown  in  Fig.4.18. 

Our  active  polygon  technique  is  also  suitable  for  document  image  segmentation  because  a  doc¬ 
ument  page  usually  has  two  different  kinds  of  texture,  namely  text  and  images,  which  have  distinct 
probability  distributions.  We  show  two  examples  of  document  segmentation:  Fig.4.19  with  a  single 
active  polygon,  and  Fig.4.20  with  two  active  polygons. 

4.6.3  Conclusions 

In  this  chapter,  we  have  presented  a  polygon  propagation  model  to  capture  particularly  textured 
objects  in  images.  A  new  ODE  model  was  developed  to  move  polygon  vertices.  In  addition,  a 
new  global  polygon  regularizer  was  introduced  to  avoid  degeneracy  during  polygon  propagation. 
Adaptation  of  a  favorable  divergence  measure,  the  Jensen-Shannon  divergence,  as  an  integral  form 
(energy  functional)  of  our  ODEs  lead  to  quite  a  powerful  unsupervised  texture  segmentation  tech¬ 
nique  which  was  validated  by  numerical  results.  The  only  assumption,  which  is  valid  in  many 
applications,  in  our  model  on  the  target  image  regions  is,  that  they  should  be  simply  connected. 
This  is  on  account  of  most  of  the  cases,  whether  in  natural  images  as  in  zebra,  or  man-made  object 
images  as  in  airplanes,  texture  regions  are  simply  connected. 
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Figure  4.17:  Each  airplane  is  captured  by  a  handful  of  vertices. 
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Figure  4. 18:  A  submarine  figure  is  captured  by  a  handful  of  vertices. 
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Figure  4.19:  An  active  polygon  nicely  segments  a  document  image  scanned  from  an  article. 
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Figure  4.20:  Two  active  polygons  capture  text  and  image  regions  of  a  page  from  an  article. 
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Chapter  5 


Applications  of  Active  Polygons 


Often,  the  goal  of  low  level  vision  is  to  extract  significant  features,  e.g.  boundaries  of  objects, 
from  an  image,  for  a  subsequent  use  in  recognition  by  the  high  level  vision  stage.  In  this  chapter, 
we  present  applications  tailored  to  the  active  polygon  model  that  we  have  developed  for  extracting 
compact  descriptions  of  objects  in  an  image.  This  representation  is  shown  to  be  useful  in  higher 
level  applications.  With  a  resulting  handful  of  vertices  for  an  object  description,  our  method  is 
highly  flexible,  and  may  be  efficiently  applied  to  object  tracking  in  video  sequences,  and  object 
recognition  in  visual  information  retrieval  context: 

1.  Video  Object  Tracking:  We  present  a  technique  for  object  tracking  in  image  sequences, 
which  makes  the  use  of  the  active  polygon  framework  developed  in  Chapter  4.  In  this  ap¬ 
proach,  upon  capturing  boundaries  of  an  object  by  a  polygon  from  the  first  few  frames  of 
the  image  sequence,  motion  segmentation  of  a  polygonal  object  can  be  quickly  achieved  by 
invoking  only  the  vertex  locations  and  the  adjacent  edges  in  the  formulation.  We  carry  out  a 
velocity  field  estimation  at  an  active  polygon  vertex  using  the  optical  flow  constraint  on  its 
two  adjacent  edges.  A  spatial  segmentation  phase  follows  to  further  refine  the  object’s  vertex 
locations  estimated  by  the  optical  flow.  The  advantage  of  our  region-based  active  polygons 
over  continuous  active  contours  in  object  tracking  in  video  applications  is  highlighted  by  its 
provision  of  a  compact  representation  of  object  features,  particularly  for  simply  connected 
target  shapes,  hence  will  be  essential  for  their  tracking  [143]. 

2.  Object  Recognition:  We  present  a  technique  for  object  recognition  which  makes  use  of  the 
compact  shape  representation  extracted  from  an  image  by  the  active  polygons.  A  signature 
function  computed  from  the  polygonal  representation  is  matched  with  the  signatures  of  the 


model  objects  in  a  small  library,  via  a  cross-correlation  operation.  An  intra-class  retrieval 
experiment,  using  photographs  of  model  airplanes  is  carried  out,  and  it  is  shown  that  active 
polygons  can  be  successfully  applied  to  such  problems. 

These  two  applications  are  explored  in  the  remainder  of  this  Chapter. 


5.1  Active  Polygons  for  Object  Tracking 

5.1.1  Overview  on  Video  Object  Tracking  Methods 

Recent  developments  in  digital  technology  have  increased  acquisition  of  digital  video  data,  which 
in  turn  have  led  to  more  applications  in  video  processing.  Video  sequences  provide  additional 
information  about  how  scenes  and  objects  change  over  time  when  compared  to  still  images.  The 
problem  of  tracking  moving  objects  remains  of  great  research  interest  in  computer  vision  on  account 
of  various  applications  in  video  surveillance,  monitoring,  robotics,  and  video  coding.  For  instance, 
MPEG-4  video  standard  introduced  video  object  plane  concept,  and  a  decomposition  of  sequences 
into  object  planes  with  different  motion  parameters  [73].  Video  surveillance  systems  are  needed  in 
traffic  and  highway  monitoring,  in  law  enforcement  and  security  applications  by  banks,  stores,  and 
parking  lots.  Algorithms  for  extracting  and  tracking  over  time  moving  objects  in  a  video  sequence 
are  hence  of  importance. 

Tracking  methods  may  be  classified  info  fwo  cafegories  [114]:  (i)  Motion-based  approaches, 
which  use  motion  segmentation  of  temporal  image  sequences  by  grouping  moving  regions  over 
time,  and  by  estimating  their  motion  models  [10,19,99,145].  This  region  tracking,  not  being 
object-based  is  not  well-adapted  to  the  cases  where  a  prior  shape  knowledge  of  the  moving  ob¬ 
ject  is  provided,  (ii)  Model-based  approaches  track  objects  using  a  template  of  the  3D  object  such 
as  3D  models  in  [50,79,91,95, 120].  Usage  of  this  high  level  semantic  information,  yields  robust 
algorithms  at  a  high  computational  complexity  cost. 

Another  classification  of  object  tracking  methods  due  to  [114]  is  based  on  the  type  of  infor¬ 
mation  that  the  tracking  algorithm  uses.  Along  these  lines,  tracking  methods  which  exploit  either 
boundary-based  information  or  region-based  information,  have  been  proposed: 

1.  Boundary-based  methods  use  the  boundary  information  along  the  object’s  contour,  and  are 
flexible  because  usually  no  objecf  shape  model  and  no  mofion  model  are  required.  However, 
methods  using  snake  models  such  as  [12,49,67],  employ  parameterized  snakes  (such  as  B- 
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splines),  and  constrain  the  motion  by  assuming  certain  motion  models,  e.g.  rigid,  or  affine. 
In  [12],  a  contour’s  placement  in  a  subsequent  frame  is  predicted  by  an  iterative  registration 
process,  where  rigid  objects,  and  rigid  motion  are  assumed.  In  another  tracking  method  with 
snakes  [89],  the  motion  estimation  step  is  skipped,  and  the  snake  position  from  any  given 
image  frame  is  carried  to  the  next  frame.  Other  methods  employ  geodesic  active  contour 
models  [24]  (which  also  assumes  rigid  motion  and  rigid  objects),  and  [114]. 

2.  Region-based  methods  such  as  [10,99, 145]  segment  a  temporal  image  sequence  into  re¬ 
gions  with  different  motions.  Regions  segmented  from  each  frame  by  a  motion  segmentation 
technique  are  matched  to  estimate  motion  parameters  [9].  They  usually  employ  parametric 
motion  models,  and  they  are  computationally  more  demanding  than  boundary -based  tracking 
methods  because  of  the  cost  of  matching  regions. 

Another  tracking  method,  referred  to  as  Geodesic  Active  Regions  [113],  incorporates  both 
boundary-based  and  region-based  approaches.  An  affine  motion  model  is  assumed  in  fhis  fechnique, 
and  ifs  compufafional  cosf  is  significanf  on  accounf  of  many  differenf  esfimafion  sfeps  involved. 

Our  goal  is  fo  build  on  fhe  existing  achievemenfs,  and  fhe  corresponding  insighf  fo  develop  a 
simple  and  efficienf  boundary-based  fracking  algorifhm  well  adapfed  fo  polygonal  objecfs.  This  is 
in  effecf  an  exfension  of  our  evolution  models  which  use  region-based  dafa  disfribufions  fo  capfure 
polygonal  objecf  boundaries  . 

5.1.2  Motion  Estimation 

Mofion  of  objecfs  in  3D  real  world  scene  are  projecfed  onfo  2D  image  plane,  and  fhis  projecfed 
mofion  is  referred  fo  as  “apparenf  mofion”,  or  “2D  image  mofion”,  or  sometimes  as  “optical  flow”, 
which  is  fo  be  esfimafed.  In  a  time-varying  image  sequence,  I{x,y,t)  :  [0,  a]  x  [0,  6]  x  [0,  T]  — )■  M+ , 
image  mofion  may  be  described  by  a  2-D  vector  field  V{x,  y,t),  which  specifies  fhe  direction  and 
speed  of  fhe  moving  fargef  af  each  poinf  {x,y).  The  measuremenf  of  visual  mofion  is  equivalenf 
fo  computing  V{x,y,t)  from  I{x,y,t)  [60].  Estimating  fhe  velocify  field  remains  an  imporfanf 
research  fopic  in  lighf  of  ifs  ubiquitous  presence  in  many  applications  and  as  rellecfed  by  fhe  wealfh 
of  previously  proposed  fechniques. 

The  mosf  popular  group  of  mofion  esfimafion  fechniques  are  referred  fo  as  differential  tech¬ 
niques,  and  solve  an  optical  flow  equafion  which  slates  lhal  intensify  or  brighlness  of  an  image 
remains  consfanf  wilh  lime.  They  use  spalial  and  temporal  derivalives  of  fhe  image  sequence  in  a 
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gradient  search,  and  sometimes  referred  to  as  gradient-based  techniques.  The  basic  assumption 
that  a  point  in  the  3D  shape,  when  projected  onto  the  2D  image  plane,  has  a  constant  intensity  over 
time,  may  be  formulated  as  (x  =  {x,  y)), 

I{x  ,t)  K,  I (x  +  5x  +  5t) 

where  5x  is  the  displacement  of  the  local  image  region  at  {x  ,  t)  after  time  5t.  A  first-order  Taylor 
series  expansion  on  the  right-hand-side  yields 

I{x  ,t)  =  I(x  ,t)  +  V/  •  5x  +  IfSt  + 

where  denotes  second  and  higher  order  derivatives.  Dividing  both  sides  of  the  equation  by  5t, 
and  neglecting  O^,  the  optical  flow  constraint  equation  is  obtained  as 

^  +  VI-V=0.  (5.1) 

This  constraint  is,  however,  not  sufficient  to  solve  for  both  components  of  V  {x  ,t)  =  {u{x  ,t),v{x  ,t)), 
and  additional  constraints  on  the  velocity  field  are  required  fo  address  fhe  ill-posed  nafure  of  fhe 
problem. 

The  direction  of  mofion  of  an  objecf  boundary  B  monitored  fhrough  a  small  aperfure  A  (small 
wifh  respecf  fo  fhe  moving  unif)  (see  Figure  5.1)  can  nol  be  defermined  uniquely  (known  as  fhe 
aperture  problem).  Experimentally,  it  can  be  observed  that  when  viewing  the  moving  edge  B 
through  aperture  A,  it  is  not  possible  to  determine  whether  the  edge  has  moved  towards  the  direction 
c  or  direction  d.  The  observation  of  the  moving  edge  only  allows  for  the  detection  and  hence 
computation  of  the  velocity  component  normal  to  the  edge  (vector  towards  n  in  Figure  5.1),  with 
the  tangential  component  remaining  undetectable.  Uniquely  determining  the  velocity  field  hence 
requires  more  than  a  single  measurement,  and  it  necessitates  a  combination  stage  using  the  local 
measurements  [139].  This  in  turn  means  that  computing  the  velocity  field  involves  regularizing 
constraints  such  as  its  smoothness  and  other  variants. 

Horn  and  Schunck,  in  their  pioneering  work  [63],  combined  the  optical  flow  constraint  with  a 
global  smoothness  constraint  on  the  velocity  field  to  define  an  energy  functional  whose  minimiza¬ 
tion 

argmin  /  [(VI  •  V  +  kf  +  XW\Vu\f  +  \\Vv\\‘^)]dx  , 

Jn 
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Figure  5.1:  The  aperture  problem:  when  viewing  the  moving  edge  B  through  aperture  A,  it  is  not 
possible  to  determine  whether  the  edge  has  moved  towards  the  direetion  c  or  direetion  d. 

ean  be  earried  out  by  solving  its  gradient  deseent  equations.  A  variation  on  this  theme,  would  adopt 
an  Li  norm  smoothness  eonstraint,  (in  eontrast  to  Horn-Sehunek’s  L2  norm),  on  the  veloeity  eom- 
ponents,  and  was  given  in  [82].  Lueas  and  Kanade,  in  eontrast  to  Horn  and  Sehunek’s  regularization 
based  on  post-smoothing,  minimized  a  pre-smoothed  optieal  eonstraint 

[  W\x)[VI{x,t)-V  +  It{x,t)?, 

Jr 

where  W  {x  )  denotes  a  window  funetion  that  gives  more  weight  to  eonstraints  near  the  eenter  of 
the  neighborhood  R  [92]. 

Imposing  the  regularizing  smoothness  eonstraint  on  the  veloeity  over  the  whole  image  leads  to 
oversmoothed  motion  estimates  at  the  diseontinuity  regions  sueh  as  oeelusion  boundaries  and  edges. 
Attempts  to  reduee  the  smoothing  effeets  along  steep  edge  gradients  ineluded  modiheations  sueh  as 
ineorporation  of  an  oriented  smoothness  eonstraint  by  [105],  or  a  direetional  smoothness  eonstraint 
in  a  multiresolution  framework  by  [47].  Hildreth  [60]  proposed  imposing  the  smoothness  eonstraint 
on  the  veloeity  field  only  along  eontours  extraeted  from  time- varying  images.  One  advantage  of 
imposing  smoothness  eonstraint  on  the  veloeity  field  is  fhaf  if  allows  for  fhe  analysis  of  general 
elasses  of  motion,  i.e.,  if  ean  aeeounf  for  fhe  projeefed  motion  of  3D  objeefs  fhaf  move  freely  in 
spaee,  and  deform  over  time  [60]. 

Spatio-temporal  energy-based  methods  make  use  of  energy  eoneenfrafion  in  3D  spatio-temporal 
frequeney  domain.  A  translating  2D  image  pattern  transformed  to  the  Fourier  domain  shows  that 
its  veloeity  is  a  funetion  of  its  spatio-temporal  frequeney  [11].  A  family  of  Gabor  filters  whieh 
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simultaneously  provide  spatio-temporal  and  frequeney  loealization,  were  used  to  estimate  veloeity 
eomponents  from  the  image  sequenees  [46, 59]. 

Correlation-based  methods  estimate  motion  by  eorrelating  or  by  matehing  features  sueh  as 
edges,  or  bloeks  of  pixels  between  two  eonseeutive  frames  [135],  either  as  bloek  matehing  in  spatial 
domain,  or  phase  eorrelation  in  the  frequeney  domain.  Similarly,  in  another  elassifieation  of  motion 
estimation  teehniques,  token-matching  schemes,  first  identify  features  sueh  as  edges,  lines,  blobs 
or  regions,  and  then  measure  motion  by  matehing  these  features  over  time,  and  deteeting  their 
ehanging  positions  [139].  There  are  also  model-based  approaches  to  motion  estimation,  and  they 
use  eertain  motion  models.  Mueh  work  has  been  done  in  motion  estimation,  and  the  interested 
reader  is  referred  to  [11, 128, 133, 135]  for  a  more  eompulsive  literature. 

5.1.3  Our  Approach 

The  key  idea  in  our  proposed  teehnique  for  video  sequenees  is  eentered  around  traeking  a  relatively 
few  vertiees  together  with  their  eorresponding  edges,  whieh  in  turn  yields  a  bookkeeping  simplie- 
ity  and  henee  effieieney.  For  a  speeifieally  fast  objeet  traeking  goal,  it  may  be  suffieient  to  only 
measure  eertain  properties  of  the  veetor  field  V  (x,  y,  t),  as  a  more  eomplete  and  preeise  eharae- 
terization  is  rather  redundant.  Moreover,  a  very  sparse  set  of  motion  measurements  only  along  the 
eontour  points,  is  required  in  lieu  of  a  dense  veloeity  field  estimation  over  the  whole  image.  A 
veloeity  estimation  is  proposed  in  tandem  with  a  fast  spatial  and  motion  segmentation  whieh  fol¬ 
low  the  initial  objeet  delineation  in  the  first  few  frames.  We  earry  out  a  veloeity  field  estimation 
at  an  aetive  polygon  vertex  using  the  optieal  flow  eonstraint  on  its  two  adjaeent  edges.  A  spatial 
segmentation  phase  follows  for  a  further  refinement  of  the  objeet’s  vertex  loeations  estimated  by  the 
optieal  flow.  As  diseussed  below,  this  approaeh  presents  several  advantages  over  aetive  eontours  in 
video  traeking;  not  the  least  of  whieh  is  its  feature  parsimony  whieh  in  turn,  greatly  simplifies  any 
traeking  applieation. 

Inspired  by  our  spatial  image  segmentation  on  a  still  image,  we  aehieve  an  ineorporation  of  a 
veloeity  field  estimation  in  ease  of  an  image  sequenee  as  deseribed  next. 
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5.1.4  Velocity  Estimation  At  Vertices 

Recall  that  our  polygon  propagating  ODEs  written  for  each  vertex  P  ^  laid  on  a  single  image  are 
given  by 

^  Pf{Lip,Pk-i,P  k))dp 

+  [  i'^-p)f{Lip,Pk,Pk+i))dp,  (5.2) 

Jo 

where  u^f.  (resp.  denotes  the  outward  unit  normal  of  edge  {Pk-i  —  Pk)  (resp.  [Pk  — 

Pk+i)),  and  t'  denotes  the  evolution  time  for  the  differential  equation.  Our  polygonal  flows  essen¬ 
tially  integrate  spatial  image  information  along  two  adjacent  edges  of  a  vertex  Pk  to  determine  its 
speed  and  direction.  This  is  exactly  the  same  idea  we  will  use  in  the  estimation  of  a  velocity  field  at 
a  vertex  of  the  active  polygon  laid  on  a  time-varying  image  sequence.  Our  goal  is  to  estimate  the  ve¬ 
locity  vector  at  each  vertex  Pk  using  the  two  adjacent  edges  as  shown  in  Fig.  5.2.  The  velocity  field 


T 


k-1 

Figure  5.2:  2-D  velocity  field  along  two  neighbor  edges  of  a  polygon  vertex. 

V  {x,y)  at  each  point  on  an  edge  may  be  represented  as  V  (p)  =  v-^{p)u-^{p)  +  v'^{p)u'^{p), 
where  u'^{p)  and  u-^{p)  are  unit  vectors  in  the  directions  tangent  and  normal  to  an  edge.  The 
optical  flow  constraint  given  by 

+  V/  •  V  =  0,  (5.3) 

ot 

provides  a  way  to  estimate  the  component  of  the  velocity  field  directly  from  the  time-varying 
image  I{x,y,t).  Once  an  active  polygon  locks  onto  a  target  polygonal  object,  the  unit  direction 
vectors  u  and  tt  ^  are  also  known  immediately.  A  set  of  local  measurements  may  easily  be 
obtained  from  the  image  brightness  constraint 

v-^ix,y)  =  -Itix,y)/\\VI{x,y)\\,  (5.4) 
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which  gives  the  magnitude  of  velocity  field  in  the  direction  orthogonal  to  the  local  edge  structure. 
This  initial  pointwise  measurement  is  preferred,  because  different  parts  of  the  moving  object  may 
have  different  motions  under  general  conditions  such  as  unrestricted  motion  [139].  However,  instan¬ 
taneous  measurements  are  insufficient  to  determine  the  motion,  and  the  local  measurements  must 
somehow  be  integrated  for  a  better  estimation  of  the  velocity  at  a  point.  Recall  that  this  was  also 
motivated  by  the  aperture  problem  which  showed  the  necessity  of  a  combination  stage  over  several 
measurements  on  the  image.  For  a  fast  estimation  of  a  velocity  at  a  polygon  vertex,  we  utilize  a 
joint  contribution  from  two  edges  of  a  vertex  to  infer  the  resultant  motion  at  the  vertex.  Expecting 
a  sensitivity  of  such  instantaneous  normal  velocity  measurements  to  noise,  we  integrate  them  in  a 
weighted  manner  along  two  neighboring  edges  of  a  vertex  to  yield  a  more  robust  estimation.  This 
leads  us  to  write  a  velocity  estimation  scheme  at  each  vertex  of  an  active  polygon  as 


dPk 

dt 


=  Vk 


=  u 


+  U 


l,k 


fo  PV'^{L(j),P  k-i,P  k))dp 
fo  P  (^P 

fd(i  -p)  v~^(L(p,Pk,Pk+i))dp 
fd(l-p)  dp 


(5.5) 


for  A:  =  1,  ...ra.  To  introduce  further  robustness  and  to  achieve  more  reliable  estimates  in  the  course 
of  computing  w*-,  we  may  make  use  of  smoother  spatial  derivatives  (larger  neighborhoods)  [43, 62, 
82]. 

To  proceed  with  tracking  of  a  polygonal  object,  we  first  apply  the  previously  described  segmen¬ 
tation  technique  to  initially  delineate  the  moving  object  boundaries.  Upon  initialization,  we  assume 
that  an  object  of  interest  has  been  already  outlined.  Using  Eq.  (5.5),  the  velocity  vector  is  estimated 
from  two  images  I{x,y,t)  and  I{x,y,t  +  1),  and  is  subsequently  utilized  to  move  the  active  poly¬ 
gon’s  handful  of  vertices  to  new  locations  on  the  next  image  I{x,y,t  +  1)  in  the  sequence.  The 
ODE  in  Eq.  (5.2)  is  then  run  for  a  short  time  for  further  refinement  of  the  polygon  delineation  of  the 
moving  object.  We  substantiate  our  approach  by  the  tracking  examples  given  in  the  next  section. 


5.1.5  Experimental  Results 

Object  Tracking  in  IR  (Infra  Red)  Video  Sequences:  A  very  suitable  application  of  our  object 
tracking  method  is  to  target  tracking  in  IR  video  sequences.  In  an  IR  image,  a  target  usually  appears 
as  a  slightly  bright  or  dark  spot,  most  of  the  time  in  the  form  of  a  rough  polygonal  shape,  when 
compared  to  the  background  terrain.  Our  active  polygons  hence  are  perfectly  suited  to  locking  onto 
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such  targets  and  tracking  them  by  the  fast  method  we  have  introduced  in  the  previous  section. 

We  show  in  Fig.  5.3,  and  Fig.  5.4,  several  snapshots  (progression  from-left-to-right,  top-bottom 
in  all  figures)  from  two  IR  sequences  respectively,  in  which  the  targets,  a  bright  triangular  shape 
in  Fig.  5.3,  and  a  bright  rectangular  shape  in  Fig.  5.4,  are  successfully  tracked  by  our  technique. 
The  active  polygon  at  each  snapshot  is  depicted  in  yellow  (white  in  b/w  print-outs).  Similarly,  in 
Fig.  5.5,  a  target  which  is  probably  a  tank,  is  tracked  as  can  be  observed  in  the  given  snapshots,  and 
the  final  represenfafion  provided  by  fhe  active  polygon  closely  resembles  a  lank-like  slruclure. 

Here,  if  should  be  emphasized  fhal  fhe  gain  of  our  approach  is  again  fwo-fold  because  in  addifion 
lo  fracking,  fhe  shape  of  fhe  Iracked  objecl  is  also  available  wilh  a  very  compacl  description,  as  a 
bonus  fealure. 

In  fhe  nexl  IR  sequence  example  (Fig.  5.6),  fhe  Iracker  approaches  gradually  lo  a  possible  large! 
sile,  which  appears  as  a  brighl  region,  and  fhe  polygon  adapls  lo  Ihe  oulline  of  Ihis  largel  region 
during  fracking.  As  a  Iasi  example  lo  IR  sequence  fracking.  Fig.  5.7  shows  a  quite  noisy,  cluttered 
terrain,  where  again  a  largel  site  appears  as  a  dark  spotted  foreground  region,  highly  embedded  in 
Ihe  background  terrain  clutter.  The  active  polygon  finds  Ihe  roughly  Iriangular  oulline  of  Ihis  site 
during  ils  successful  fracking  Ihrough  Ihe  IR  image  sequence. 

In  Ihe  nexl  example  (Fig.  5.8),  we  Iracked  a  toy  airplane,  a  model  of  an  F22,  which  is  Iranslaled 
on  a  conveyor  bell.  Note  lhal  once  Ihe  active  polygon  locks  onto  Ihe  largel  objecl,  il  follows  Ihe 
movemenl  of  Ihe  largel  unaffected  by  Ihe  camera  movemenls  such  as  zoom,  pan,  Iranslalion,  and 
rolalion.  Similarly,  a  model  rockel  which  is  aboul  to  be  launched,  is  zoomed  oul  by  Ihe  camera  in 
Fig.  5.9,  and  is  well  Iracked  by  our  algorilhm. 

The  Iasi  Iwo  fracking  examples  show  a  potential  application  of  our  technique  to  Iraflic  moni¬ 
toring  where  overhead  cameras  can  give  information  on  slate  of  Iraffic,  ils  density,  Ihe  speed  of  Ihe 
vehicles,  or  abnormal  silualions  on  Ihe  road.  In  Ihe  firsl  sequence,  shown  in  Fig.  5.10,  a  fixed  cam¬ 
era  has  been  mounted  on  a  hilltop  overviewing  Ihe  Tryon  Road,  Raleigh,  NC.  A  black  automobile 
entering  Ihe  scene  from  Ihe  left,  is  successfully  Iracked  by  an  active  polygon  until  il  leaves  Ihe  scene 
al  Ihe  olher  end  of  Ihe  road  hence  oul  of  Ihe  camera  range.  In  anolher  road  sequence  (of  COST  211 
projecl  [1]),  shown  in  Fig.  5.11,  a  Iruck  approaching  from  Ihe  horizon  on  Ihe  road  (left),  is  Iracked 
via  an  active  polygon.  Active  polygons  are  successfully  applied  to  fracking  of  motor  vehicles,  be¬ 
cause  Ihe  shapes  of  Ihese  man-made  vehicles  are  polygonal,  making  Ihe  motion  fracking  Ihrough 
vertices  of  a  polygonal  shape  feasible  and  practical. 
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Figure  5.3:  Tracking  of  an  IR  image  sequence  I. 
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Figure  5.4:  Tracking  of  an  IR  image  sequence  II. 
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Figure  5.5:  Tracking  of  an  IR  image  sequence  III. 
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Figure  5.6:  Tracking  of  an  IR  image  sequence  IV. 
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Figure  5.7:  Tracking  of  another  IR  image  sequence  V. 
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Figure  5. 10:  Tracking  of  the  black  automobile,  entering  the  scene  in  the  top  left  figure,  and  leaving 
the  scene  in  the  bottom  right  figure  (sequence  provided  by  courfesy  of  Alper  Unal,  Civil  Engineering 
Depf,  NCSU). 


Figure  5.11:  Tracking  of  a  fruck  approaching  from  the  far  left  (sequence  of  COST211  project). 
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5.2  Object  Recognition  for  Visual  Information  Retrieval 


Visual  information  retrieval  refers  to  the  problem  of  retrieving  images  or  image  sequenees  from  a 
database,  whieh  are  relevant  to  a  query.  Visual  information  retrieval  problem  arose  as  a  new  topie 
of  researeh  on  aeeount  of  emerging  multimedia  applieations,  the  availability  of  large  image  and 
video  arehives,  and  one’s  ability  of  sharing  and  distributing  image/video  data  over  large  bandwidth 
eomputer  networks  [15].  First  generation  visual  information  retrieval  systems  allowed  querying  by 
strings  or  text  with  whieh  eapturing,  for  instanee,  pereeptual  elements  of  a  shape  or  visual  features 
sueh  as  the  outline  of  an  objeet  is  diffieult.  In  the  new  generation  visual  information  retrieval 
systems,  a  more  vital  description  of  visual  content  in  agreement  with  human  perception  is  adopted, 
where  image  processing  and  computer  vision  techniques  are  essential  for  automatic  extraction  of 
visual  features  like  color,  shape,  and  texture,  from  the  image  data. 

The  most  recent  member  of  the  MPEG  family,  MPEG-7,  specifies  standardization  of  visual 
descriptors  (e.g.  color,  shape,  texture,  and  motion),  and  description  schemes  (complex  descriptions 
specifying  structure  and  semantic  relationships  among  descriptors)  of  multimedia  information  to 
allow  efficient  and  fast  content-based  access  and  search.  Development  of  visual  retrieval  systems  is 
hence  not  only  gaining  importance  but  becoming  a  very  important  problem  as  well. 

Eor  content-based  image  retrieval,  particularly  when  objects  in  an  image  are  of  interest,  object 
recognition  is  a  desirable  goal,  where  a  decision  of  whether  or  not  an  observed  object  corresponds 
to  a  model  in  the  database  has  to  be  made.  In  the  traditional  matching  approaches,  a  comparison 
between  the  query  image  and  the  database  images  is  made,  and  an  ordering  according  to  the  mea¬ 
sured  similarity  is  output.  Selecting  a  similarity  measure  is  very  important,  and  depends  on  the 
representation  of  perceptual  features. 

One  of  the  most  meaningful  descriptions  of  an  object  is  through  its  shape,  and  a  simple  and 
efficient  way  for  representing  a  shape  is  through  a  set  of  features  extracted  from  the  image.  Eeatures 
such  as  the  area,  the  elongatedness,  the  circularity,  the  orientation  of  the  major  axis  describe  global 
shape  characteristics,  whereas  features  like  corners,  characteristic  points  on  the  boundary  describe 
local  characteristics.  Eocal  features  are  less  susceptible  to  noise  due  to  a  partial  description  they  may 
provide  even  in  the  presence  of  occlusions.  Representing  a  shape  through  features  that  describe  the 
boundary  of  the  object,  either  locally  or  globally,  is  a  popular  method.  Considering  set  of  all  pixels 
on  the  boundary  is,  however,  not  efficient  as  the  number  of  points  is  usually  too  large.  Reduction 
of  the  dimension  of  the  shape’s  boundary  representation  by  focusing  on  perceptually  salient  points 
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(e.g.  curvature  extrema,  corner  points),  is  desirable.  Approximating  shapes  through  a  polygonal 
curve  has  previously  been  mentioned  in  shape  retrieval  techniques  [98].  In  [98],  however,  shapes 
are  already  input  as  a  polygon  to  the  retrieval  system,  and  polygonal  approximation  methodology 
is  not  discussed  in  the  case  of  a  continuous  input  shape.  A  group  of  interest  points,  i.e.  vertices, 
form  features  which  are  encoded  and  collectively  form  a  feature  index  for  the  query  shape.  Our 
approach  is  taking  rather  than  an  indexing  approach,  a  shape  matching  approach  (will  be  explained 
in  the  next  subsection),  in  an  image  retrieval  problem  which  is  more  involved  than  a  shape  retrieval 
problem.  Also  in  [98],  a  discrimination  between  object  recognition  problem  and  similar  shape 
retrieval  problem  has  been  made,  and  the  former  is  defined  to  require  an  automatic  recognition 
of  an  unknown  object  (without  user  intervention),  whereas  the  latter  one  is  assumed  to  require 
user  assistance.  We  do  not  make  such  distinction  between  the  two  problems  here,  and  assume  an 
automatic  object  recognition  scenario  in  both  cases. 

Model-based  Object  Recognition  using  invariants 

Existing  object  recognition  systems,  which  usually  address  the  problem  of  selecting  the  object 
model  that  best  matches  the  observed  image,  are  referred  to  as  model-based  [3,21,84, 104].  In 
these  systems,  the  model  objects,  and  the  observed  object  are  described  by  a  set  of  feature  points, 
whose  extraction  is  the  first  task  to  be  completed.  This  is  usually  followed  by  a  comparison  where 
a  quantitative  similarity  between  two  shapes,  is  carried  out  by  a  distance  measure  between  their  de¬ 
scriptions.  The  computed  shape  descriptions  should  be  robust  in  the  context  of  invariance  property 
so  crucial  to  reliable  object  recognition. 

Geometric  invariants  are  properties  of  geometric  configurations  of  features  which  remain  un¬ 
changed  under  a  certain  class  of  transformations.  Geometric  characteristics  of  objects  are  usually 
distorted  by  a  perspective  projection  of  a  scene  from  3D  space  (real  world)  to  a  2D  image  plane.  For 
instance,  length  is  an  invariant  under  rigid  transform  by  rotation  and  translation,  but  it  is  changed 
under  a  more  general  class  of  transformations  like  affine,  or  perspective  fransformafion.  Parallelism 
is  preserved  under  an  affine  fransformafion  (e.g.  a  square  is  fransformed  fo  an  arbifrary  parallel¬ 
ogram).  On  fhe  ofher  hand,  collinearify  of  poinfs  is  preserved  under  a  projective  fransformafion 
(equivalenf  fo  perspective  viewing)  [44,57].  Imporfanf  invarianfs  of  affine  fransformafion  are  par¬ 
allel  lines,  rafio  of  lengfhs  of  parallel  line  segmenfs,  and  ratio  of  areas.  The  mosf  fundamenlal 
invarianf  of  a  projective  fransformafion  is  fhe  cross  rafio  of  four  collinear  poinfs  (a  rafio  of  ratios  of 
lengfhs  on  a  line). 
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In  model-based  vision,  affine  geometry  is  often  assumed  beeause  of  the  availability  of  a  larger 
number  of  invariant  features,  and  the  sufticieney  of  features  neeessary  to  uniquely  eompute  model 
pose  under  an  affine  projeetion  rather  than  a  perspeetive  projeetion. 

The  erueial  step  of  invariant-based  reeognition  is  in  the  reduetion  of  the  huge  amount  of  in¬ 
formation  in  an  image  to  a  small  number  of  features  whieh  remains  invariant  to  viewing  angle, 
oeelusion,  lighting  eonditions  whieh  are  in  turn  quintessential  to  robustness. 

There  are  two  approaehes  to  model-based  objeet  reeognition  using  geometrie  invariants. 

•  Search-based  matching  approach:  Invariant  features  are  eomputed  from  the  observed  im¬ 
age  and  are  eompared  with  pre-eomputed  invariant  features  of  all  models  in  the  database.  This 
exhaustive  seareh  teehnique  is  good  for  small  databases,  but  problematie  for  large  databases. 

•  Indexing  approach:  The  derived  transformation-invariant  deseription  is  used  as  a  key  in 
indexing  a  database  of  objeet  models  for  reeognition.  In  the  so-ealled  Geometrie  Hashing 
teehnique  [85],  objeets  are  represented  as  sets  of  geometrie  features  sueh  as  points  or  lines, 
and  their  geometrie  relations  are  eneoded  using  minimal  feature  sets  to  index  a  hash  table 
prepared  beforehand  in  a  learning  phase.  In  the  reeognition  phase,  a  vote  of  the  hits  to  entries 
of  the  hash  table  by  the  feature  indiees  of  the  observed  image  are  tallied,  and  a  database  model 
objeet  that  reeeives  the  maximum  number  of  hits  is  ehosen  as  the  eorresponding  objeet  model 
for  the  queried  unknown  objeet.  These  methods  are  preferred  for  effieiently  aeeessing  very 
large  model  databases. 

Another  elassitieation  in  image  retrieval  problem  can  be  made  in  terms  of  the  homogeneity  of 
the  database.  The  inter-class  retrieval  problem  involves  distinguishing  different  classes  of  objects 
like  airplanes,  animals,  and  cars.  The  intra-class  retrieval  problem  on  the  other  hand  involves 
distinguishing  objects  from  the  same  class,  for  instance  airplane  types  F16  and  FI  17.  An  intra-class 
indexing  application  to  digital  fish  databases  is  given  in  [4].  Inter-class  retrieval  may  be  handled 
by  comparing  more  global  features  of  a  shape,  while  intra-class  retrieval,  however,  requires  more 
sophisticated  approaches  to  feature  comparison. 

In  the  next  subsection,  we  present  an  intra-class  image  retrieval  study  we  have  carried  out  as  an 
application  of  our  active  polygon  model. 
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Figure  5.12:  Intra-class  Image  Retrieval  Scenario. 

5.2.1  Active  Polygons  for  Intra-Class  Image  Retrieval 

Complexity  of  an  intra-class  image  retrieval  problem  arises  from  the  fact  that  objects  in  a  database, 
e.g.  different  types  of  airplanes,  and  the  query  image  are  acquired  from  different  viewing  angles,  dif¬ 
ferent  lighting  and  background  conditions.  We  performed  a  recognition  experiment  on  real  images 
as  follows.  The  object  database  contains  pictures  of  various  airplane  models.  We  took  photographs 
of  the  hve  different  airplane  models  from  a  top  view  (shown  in  Fig.  5.12),  and  from  different  viewing 
perspectives,  and  under  different  lighting  conditions,  given  in  Fig.  5.13.  This  scenario  is  depicted 
in  Fig.  5.12,  where  a  query  airplane  image  is  input  to  the  system,  with  a  question  of  identifying  its 
membership  model.  The  test  image  (i.e.  the  query),  is  viewed  under  different  lighting  and  viewing 
conditions.  The  key  idea  is  to  extract  signature  features  from  the  test  image,  and  search  from  the 
database  for  the  object  model  that  best  matches  the  query  in  terms  of  this  signature. 

Our  basis  for  solving  this  problem  is  the  active  polygon  framework  we  have  proposed  in  Chapter 
4.  One  of  the  important  pre-processing  steps  for  object  recognition  systems  consists  of  automatically 
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Figure  5.13:  Set  of  airplane  images  under  different  eonditions. 

and  systematieally  extraeting  features,  or  the  interest  points,  from  an  image.  We  believe  that  this 
issue  is  addressed  unsatisfaetorily  at  best  in  most  of  the  previous  objeet  reeognition  work  where 
manual  or  supervised  extraetion  was  earried  out.  The  aetive  polygons  we  have  developed,  on  the 
other  hand,  ean  automatieally  extraet  important  feature  points  of  man-made  objeets  from  images.  A 
polygonal  eontour,  for  eaeh  airplane  shape  in  an  image  is  for  instanee  extraeted  by  our  model. 

Affine  Invariant  Parameterization  of  the  Polygon  Vertices 

The  extraeted  feature  points  (i.e.  the  ordered  set  of  vertex  points)  are  used  in  forming  a  signa¬ 
ture  from  eaeh  shape,  and  in  order  to  aeeount  for  some  degree  of  affine  invarianee,  we  ehoose  to 
reparameterize  the  vertex  points  by  an  affine  invarianf  parameferizafion.  An  affine  fransform  is  a 
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nonsingular  linear  transformation  followed  by  a  translation:  X  =  A  X  +  B ,  where  for  a  2D 
representation  of  shapes,  A  is  a  2  x  2  matrix,  X  and  B  are  2  x  1  vectors. 

A  curve  C  (p)  =  {X{p),Y{p)),  C  :  [a,  6]  — )■  with  parameter  p  G  [a,b]  can  be  reparameter¬ 
ized  to  a  new  parameter  s  via  a  metric  g{p)  so  that  s{p)  =  The  metric  g{p)  =  |det  < 

C  ,C  p  >  \  first  proposed  in  [7],  is  a  first  order  affine  invariant  metric  for  B  =  0,  i.e.  the  translation 
component  of  the  affine  transform  is  zero.  Note  that  the  affine  transform  on  this  metric  can  be  given 
by  ds  =  det  <  C  ,  C  p  >  dp  =  det  <  AC  ,  A  C  p  >  dp  =  detA  <  C  ,  C  p  >  dp,  and  the  metric 
is  hence  linear  under  an  affine  transform  if  B  =0. 

We  note,  using  the  divergence  theorem,  that  this  new  parameterization 

ta  =  ^l\xYp-YXp\dp 
is  indeed  the  the  area  enclosed  by  the  contour  C , 


Since  the  affine  transform  linearly  changes  the  area  with  a  constant  scale  det  A ,  a  normalization  with 
respect  to  the  total  enclosed  area  of  the  contour  is  possible.  An  initial  translation  of  the  polygon  to 
the  affine  center  of  mass,  as  a  result,  makes  this  parameterization  completely  affine  invariant,  i.e. 
invariant  to  translation,  rotation,  scale,  and  shear. 

In  [7],  a  Fourier  transform  of  the  affine-invariant  parameterized  contour  was  taken,  and  invari¬ 
ants  were  synthesized  from  the  Fourier  coefficients.  Inspired  by  this  approach,  wavelet  transform 
coefficients  of  the  affine-invariant  parameterized  contour  were  used  as  invariants  in  [3].  In  pursuit 
of  exploiting  the  parsimonious  representation  provided  by  our  active  polygons,  we  resort  to  obtain  a 
signature  that  can  be  easily  matched  for  a  small  number  of  vertices.  The  cumulative  angle  function, 
or  turning  function  [8],  measures  the  angle  between  the  counterclockwise  tangent  and  the  rr-axis  as 
a  function  of  the  arc  length  s.  It  increases  with  left-hand  turns  and  decreases  with  right-hand  turns. 
For  the  polygonal  contours  that  we  obtain  for  each  image,  the  turning  function  is  piecewise  constant, 
with  increases  or  decreases  at  the  vertex  points.  Instead  of  using  the  standard  turning  function  as 
a  signature,  we  propose  to  use,  the  turning  function  of  the  polygon  versus  the  first  order  affine  arc 
length  (i.e.measure  the  angle  of  the  counterclockwise  tangent  as  a  function  of  the  affine  arc  length). 
Although  the  resulting  signature  is  not  completely  affine  invariant,  it  will  be  observed  that  such  a 
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Figure  5.14:  An  airplane  from  two  different  views. 


parameterization  gives  better  results  in  matching  polygonal  shapes  with  different  poses. 

In  Figure  5.14,  turning  function  signatures  are  shown  for  an  airplane  from  two  different  views. 

In  light  of  the  reduced  number  of  vertices  of  our  polygonal  representation,  one  may  carry  out 
a  cross-correlation  on  the  turning  functions  of  the  two  polygons  in  order  to  account  for  the  starting 
point  variability  among  shapes.  Recall  that  cross-correlation  is  a  measure  of  similarity  between 
two  different  data  sets,  computed  by  the  sum  of  the  products  between  the  two  data  sets  at  different 
lags.  The  advantage  of  our  proposed  model  in  using  this  signature  becomes  clear  as  the  required 
number  of  lags,  i.e.,  the  dimension  of  the  cross-correlation  function,  is  very  small,  (in  fact  equals 
the  minimum  of  the  number  of  vertices  of  the  two  polygons).  For  each  lag  I,  the  normalized  cross¬ 
correlation,  so-called  the  cross-correlation  coefficient,  between  two  signature  functions,  0i  and  ©2, 
is  defined  as 

EtJ0l(ia)  -  ^{0l(fa)}][02(fa  -  1)  -  ^{02(ta)}] 
{EtJ0l(U-^{0l(ia)}F[Et.02(ia-O-^{02(ta)}F}°-^'  ^  ^ 

In  Figure  5.15,  fhe  fesf  image  is  an  FI  17  whose  polygonal  represenfafion  is  exfracfed  by  an 
acfive  polygon,  and  ifs  signafure  is  computed,  and  compared  fo  fhe  signafures  of  fhe  airplanes  in 
fhe  model  image  sef  by  way  of  fhe  cross-correlafion  coefficienf.  The  model  objecf  whose  signafure 
achieves  a  maximal  cross-correlafion  wifh  fhaf  of  fhe  fesf  objecf’s  is  picked  as  a  fop  mafch.  The 
consecufive  mafches  are  also  ordered  wifh  respecf  fo  fhis  criferion  (Figure  5.16).  If  can  be  observed 
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Figure  5.15:  A  query  (an  FI  17)  to  the  model  image  set,  and  the  best  match  are  shown. 

that  the  top  match  is  indeed  the  correct  model  object  from  the  image  set. 

Another  query  image,  an  F16,  and  the  corresponding  retrieved  object  models  from  the  data  set 
is  shown  in  Figure  5.17  and  Figure  5.18.  The  best  matching  model  is  an  F16. 

A  Statistical  Approach:  We  also  carry  out  a  statistical  study  when  generating  the  signatures 
before  matching.  Due  to  an  inherent  variability  in  real  images,  the  extracted  feature  set  may  slightly 
vary  for  different  initial  conditions  of  an  active  polygon.  To  account  for  this  variability,  we  generated 
random  active  polygon  initializations  on  an  image,  and  hence  obtained  several  realizations  of  the 
target  shape  in  terms  of  its  feature  points  (see  eight  different  realizations  of  an  F16  airplane  in 
Figure  5.19).  An  average  signature  over  the  signatures  of  these  realizations  may  then  be  obtained, 
after  aligning  the  shapes  by  cross-correlating  their  turning  functions  vs.  standard  arc  length  which 
is  opted  to  be  used  in  this  experiment.  The  average  turning  function  of  the  eight  realizations  in 
Figure  5. 19,  computed  in  this  way,  is  shown  in  Figure  5.20. 

Two  retrieval  examples  using  the  match  over  averaged  turning  functions  are  shown  in  Fig¬ 
ure  5.21.  The  correct  model  from  the  data  set  is  once  again  picked  in  both  of  the  cases. 

We  note  that  the  application  of  active  polygons  to  an  intra-class  image  retrieval  problem  pre¬ 
sented  above,  requires  a  more  thorough  performance  evaluation,  with  a  larger  database  and  extensive 
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Figure  5.16:  The  next-to-top  matehes  are  shown  for  the  test  airplane  in  Fig.  5.15. 
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Figure  5.17:  A  query,  which  is  an  F16,  to  model  image  set,  and  a  top  match  are  shown. 
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Figure  5.18:  The  next-to-top  matches  are  shown  for  the  test  airplane  in  Fig.  5.17. 


Oupul  bTi^viim  ■vhhi  fdm 


Figure  5.19:  Outputs  of  several  random  initializations  for  F16  airplane. 
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Figure  5.20:  Average  turning  function  versus  arc  length  over  eight  different  realizations  shown  in 
Fig.5.19. 

simulations.  This  study  may  be  considered  as  part  of  a  future  research  direction  motivated  by  this 
thesis. 


5.3  Conclusions 

In  this  chapter,  we  aimed  at  demonstrating  the  flexibility  and  potential  uses  of  the  active  polygon 
framework,  developed  in  Chapter  4,  in  a  varied  set  of  applications  in  computer  vision.  First,  we 
presented  an  object  tracking  application,  which  provided  a  fast  way  of  outlining  a  moving  object 
through  time  in  a  video  sequence.  An  additional  gain  provided  by  our  technique  is  the  immediate 
availability  of  the  compact  shape  description  of  the  tracked  object.  We  also  presented  an  experi¬ 
mental  study  on  applying  active  polygons  to  object  recognition,  which  also  clearly  illustrated  the 
potential  of  our  technique  in  such  applications. 
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Figure  5.21:  A  query,  (a)  F16,  (b)  F22,  to  model  image  set,  and  a  top  match  using  averaged  turning 
functions  are  shown. 
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Chapter  6 


Contributions  and  Future  Research 


In  this  chapter,  contributions  of  this  thesis,  and  the  conclusions  drawn  from  the  associated  research 
work  and  the  results  are  presented.  Suggestions  for  future  research  directions  related  to  this  thesis 
are  provided. 

6.1  Contributions  of  the  Thesis 

In  this  thesis,  we  have  presented  contributions  in  the  development  of  new  image  processing  algo¬ 
rithms  based  on  curve  and  polygon  evolution  models.  In  Chapter  1,  the  use  of  curve  evolution 
techniques  for  image  processing  are  motivated,  and  then  in  Chapter  2,  the  background  material 
fundamental  to  the  curve  evolution  concept  is  given  as  preliminaries. 

6.1.1  Contributions  to  Curve  Evolutions  for  Nonlinear  Filtering 

In  Chapter  3,  we  described  a  new  approach  to  curve  and  image  smoothing  through  a  new  class  of 
curve  evolution  equations  designed  to  preserve  prescribed  polygonal  structures  in  an  image  while 
removing  noise.  These  curve  evolutions  which  are  obtained  via  a  directional  generalization  of 
the  geometric  heat  equation  that  circularizes  any  closed  curve,  are  applied  to  smooth  noisy  curves 
without  destroying  their  significant  features  on  which  a  prior  knowledge  is  assumed.  A  contribution 
in  conjunction  with  these  flows  is  a  local  stochastic  formulation  of  the  geometric  heat  equation 
providing  an  alternative  microscopic/macroscopic  view  of  this  equation,  whose  insight  led  to  pre¬ 
set  and  stabilized  polygonal  approximations.  The  theoretical  framework  of  the  new  class  of  curve 
evolutions  which  leads  to  a  nonlinear  filtering  along  known  salient  structural  directions  in  an  image 
constitutes  the  first  contribution  of  this  thesis.  Additional  contributions  come  from  the  availability 
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of  various  applications  besides  image  smoothing.  We  demonstrated  that  the  designed  flows  have  the 
ability  to  morph  an  input  shape  into  various  other  shapes,  and  thus  constructed  continuous  shape 
transformations  have  their  potential  use  in  shape  analysis  and  shape  recognition. 

6.1.2  Contributions  to  Active  Contours  for  Unsupervised  Texture  Segmentation 

Chapter  4  introduced  a  novel  variational  framework  for  an  unsupervised  texture  segmentation  prob¬ 
lem.  The  first  contribution  in  this  framework,  is  the  new  ordinary  differential  equation  model  that 
moves  polygon  vertices,  obtained  as  the  gradient  flow  of  a  generic  region-based  energy  functional 
for  an  active  contour.  Under  the  new  evolution  model,  the  evolving  contour  is  a  polygon  whose 
every  vertex  is  propagated  by  an  overall  speed  function  integrated  along  its  two  adjacent  polygon 
edges.  This  is  a  significant  improvement  over  other  snakes,  geodesic  active  contours,  and  region- 
based  active  contour  techniques.  The  latter  group  of  techniques  which  are  closer  to  our  model,  fuse 
global  information  inside  and  outside  the  curve  with  a  pointwise  local  measurement,  and  is  hence 
not  amenable  to  speed  functions  that  try  to  capture  higher-order  statistical  features  present  in  most 
of  the  textural  variations  in  an  image.  The  second  contribution  in  this  set-up  is  the  definition  of  an 
information-theoretic  measure,  based  on  the  Jensen-Shannon  divergence,  which  unfolds  the  textu¬ 
ral  information  through  the  underlying  probability  distributions  of  the  data,  moreover  using  easily 
computed  statistics.  Accounting  for  underlying  distributions  exactly  in  order  to  separate  two  differ¬ 
ent  textural  regions  in  an  image,  although  desirable,  would  involve  estimation  of  densities  possibly 
through  histograms,  would  hence  be  demanding  in  terms  of  computations.  Utilizing  an  approxi¬ 
mation  of  Shannon  entropy  in  the  definition  of  the  Jensen-Shannon  divergence,  and  in  turn  in  the 
cost  functional,  lead  to  a  fast  numerical  scheme  to  estimate  higher-order  characteristics  of  a  proba¬ 
bility  distribution  with  minimal  effort.  In  addition,  another  advantage  of  adapting  Jensen-Shannon 
divergence  as  the  energy  functional  of  active  contours  is  reflected  by  its  ability  of  generalization 
to  propagating  multiple  active  contours  whose  gradient  descent  equations  are  derived  from  a  single 
cost  functional,  hence  leading  to  a  new  coupled  set  of  active  contour  evolution  equations.  A  third 
contribution  in  the  same  active  polygon  set-up  is  the  new  polygon  regularizer,  which,  in  contrast 
to  other  polygon  evolvers,  is  global  and  uses  a  physical  basis  depending  on  electrostatics.  This  is 
also  a  significantly  new  contribution  to  the  field  of  computer  vision  in  terms  of  curve  evolutions. 
These  three  ideas  of  Chapter  4  together  result  in  a  new  active  contour  model,  which  is  significantly 
different  from  other  snake  and  active  contour  methods,  and  which  is  successfully  applied  to  more 
general  scenarios  of  images  in  terms  of  texture  and  to  various  computer  vision  applications. 
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6.1.3  Contributions  to  Computer  Vision 


Chapter  5  extended  the  eontribution  of  Chapter  4  to  time-varying  image  sequenees  for  an  objeet 
traeking  applieation.  The  additional  eontribution  here  is  an  effieient  and  simple  traeking  method¬ 
ology  whieh  builds  on  the  idea  of  traeking  a  relatively  small  number  of  vertiees  by  eomputation 
of  a  veloeity  field  at  a  vertex  along  its  two  adjaeent  edges.  We  demonstrate  through  substantiating 
examples,  that  the  aetive  polygon  framework  is  easily  adapted  to  this  seenario.  The  effieieney  of  the 
proposed  traeking  method  eomes  from  the  suffieieney  of  a  very  sparse  set  of  motion  measurements. 
The  motion  only  along  the  polygon  edges  is  used,  and  a  eoarse  estimate  at  a  vertex  integrated  in  a 
weighted  manner  the  immediately  available  loeal  measurements,  i.e.  the  eomponent  of  the  veloeity 
orthogonal  to  an  edge.  An  applieation  to  traeking  targets  in  IR  image  sequenees  validates  the  use¬ 
fulness  of  the  teehnique.  The  flexibility  of  the  aetive  polygon  framework  is  also  demonstrated  in  its 
sueeessful  applieation  to  an  experimental  objeet  reeognition  seenario. 

6.2  Future  Research  Directions 

Several  interesting  researeh  direetions  motivated  by  this  thesis  are  diseussed  next. 

6.2.1  Shape  Prior  on  the  Curve  Evolution  Model 

In  Chapter  3,  the  flows  we  designed  ineorporated  a  prior  knowledge  on  the  salient  direetions  of 
polygonal  shapes  whose  preservation  during  smoothing  is  desired.  The  form  of  the  funetions,  h{-) 
used  in  the  flows  were  trigonometrie,  and  were  based  on  the  unit  normal  angle,  i.e.  the  loeal  ori¬ 
entation  of  a  point  on  the  eurve.  They  lead  to  simple  n-gone  shapes  sueh  as  triangles,  reetangles, 
hexagons,  and  so  on,  on  aeeount  of  the  shape  prior  being  used.  An  open  question  is  how  to  ineorpo- 
rate  other  and  more  general  forms  of  shape  priors,  for  instanee  global  eharaeteristies  of  a  shape  sueh 
as  its  size,  elongatedness,  number  of  eorners  and  so  on.  A  related  question  is  that  through  what  kind 
of  loeal  and  global  features,  other  than  orientation,  more  general  shape  priors  may  be  integrated  by 
way  of  the  funetions  utilized  in  this  model.  We  allude  to  sueh  a  direetion  in  [140, 142],  where  a 
variant  of  the  h{-)  funetions  presented  in  this  thesis,  was  used  on  the  basis  of  derivative  of  the  loeal 
orientation. 
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6.2.2  A  New  Application  for  the  Curve  Evolution  Model 

In  Chapter  3,  we  elaborated  on  the  applieation  of  the  new  elass  of  eurve  evolutions  to  image  and 
shape  smoothing.  Another  interesting  applieation,  henee  an  interesting  future  researeh  direetion, 
for  these  flows  when  eonsidered  as  a  eontinuous  transformation  of  one  shape  into  another,  is  in 
eomputer  vision  and  graphies.  Morphing  in  eomputer  graphies  transforms  the  shape  of  graphieal 
objeets  sueh  as  2D  eurves,  regions,  images,  surfaees,  and  volumes.  This  may  also  be  eonsidered  as 
a  powerful  method  for  shape  analysis,  where  the  transformations  required  to  deform  one  shape  into 
the  other  provide  a  way  to  elassify  the  members  of  a  given  shape  family. 

6.2.3  A  New  Cost  Functional  for  the  Polygon  Evolution  Model 

In  Chapter  4,  the  polygon  evolution  models  we  have  developed,  form  a  new  aetive  eontour  frame¬ 
work  for  texture  segmentation  problem.  These  models  made  use  of  a  fast  numerieal  estimation 
seheme  in  eomputation  of  the  Jensen-Shannon  divergenee,  derived  from  a  first-order  approximation 
of  a  probability  density  and  the  Shannon  entropy.  An  avenue  of  future  researeh  would  be  to  de¬ 
fine  eost  funetionals  whieh  better  approximate  the  underlying  probability  distributions  of  regions, 
or  develop  more  effieient  nonparametrie  estimation  teehniques  of  sueh  measures.  Adapting  the 
same  divergenee  measure  with  a  Renyi  entropy  [56]  of  regions,  as  a  more  general  measure  than  the 
Shannon  entropy,  would  be  an  interesting  investigation  into  understanding  the  gains  that  would  be 
aehieved  in  this  ease. 

6.2.4  Constraint  on  the  Number  of  Vertices  and  a  Shape  Prior  on  the  Polygon  Evo¬ 
lution  Model 

Polygon  evolution  models  developed  in  Chapter  4,  adaptively  (by  periodieally  adding  and  removing 
vertiees)  eonform  to  an  adequate  number  of  vertiees  that  “minimally”  define  a  target  shape  in  an 
image.  The  problem  of  adding  a  eonstraint  on  the  ehoiee  of  the  number  of  vertiees  of  the  polygon 
to  the  eost  funetional  of  the  aetive  eontour,  would  be  an  interesting  future  researeh  avenue.  The 
question  of  putting  a  penalty  on  the  eost  of  adding  more  and  more  vertiees,  or  a  prior  distribution 
on  the  type  of  the  polygon,  its  shape,  and  the  number  of  vertiees,  and  finally  eoding  this  prior 
distribution  by  a  minimum  deseription  length  eriterion  would  be  a  researeh  direetion  that  holds 
promise. 


119 


6.2.5  Topological  Changes 

The  active  polygon  model  developed  in  this  thesis  makes  the  assumption  of  simply  connected  re¬ 
gions  for  the  image  segmentation  problem.  This  assumption  is  not  very  restrictive,  since  in  most  of 
the  cases,  whether  in  natural  images  as  in  zebra,  or  man-made  object  images  as  in  airplanes  (shown 
in  Chapter  4),  texture  regions  are  simply  connected.  An  investigation  in  how  to  handle  topological 
changes,  i.e.  split  and  merge  of  a  polygon,  would  be  an  improvement  to  capture  multiply  connected 
target  regions  in  an  image  for  an  extension  to  more  general  scenarios. 

6.2.6  Further  Extensions  and  Applications  for  Active  Polygons 

•  Active  polygon  methodology  is  extended  to  an  object  tracking  problem  in  Chapter  5.1.  An 
improved  motion  estimation,  by  using  more  than  the  measured  orthogonal  components  of  the 
velocity  along  two  adjacent  edges  for  each  vertex,  by  imposing  for  instance  a  smoothness 
constraint  on  the  velocity  field  along  the  contour.  This  would  enhance  the  object  tracking 
methodology,  possibly  at  a  cost  of  slowing  down  the  performance. 

•  Designing  geometric  invariants,  fully  affine  invariant  or  projective  invariant,  suitable  for  the 
output  of  the  active  polygon  model  would  result  in  a  full  extension  of  the  object  recognition 
application  framework  started  in  Chapter  5.2. 
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Appendix  A 


Stochastic  Formulation  of  a  Geometric 
Heat  Equation 


Let  us  denote  by  x  )  the  solution  to  Eq  (3.6): 


If  we  write  u{t,x)  as 

u{t,  X  )  =  X  )  +  eu{t,  X  ) 


then  eu{t,x)  satisfies 
d€.U 

—  =  fi(u2  +  eu^,u^  +  euy,u2^  +  eu^^)- 

-  {f2iu^  +  +  euy,u^y  +  eu^y)  -  /2«,u”,<j^)} 

+  hiUx+^Ux,Uy  +eUy,Uyy+eUyy)  -  f3iU^,Uy,Uyy).  ( A .  l  ) 
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Assuming  fi{-,  f2{-,  ■,  ■)  and  /3(-,  •,  •)  are  differentiable  in  their  arguments,  we  ean  expand 

in  Taylor  series  about  /2(-,-,-)  about  (u^,u^,u^y),  and  hi-rr)  about 

{u^ ,  Uy ,  Uyy ) .  For  uotattonal  simplieity,  let  us  denote  by  a: i  =  (ux ,  Uy ,  Uxx )  and  =  {u^,Uy,  u^x ) ’ 
X2  =  {ux,Uy,Uxy)  wd  x^  =  and  0:3  =  {ux,Uy,Uyy)  wd  x'^  =  {u^,u^,u^y). 

If  we  assume  that  eu{t,x)  is  small  enough,  we  ean  negleet  higher  order  terms  and  write  a  linear 
approximation  as 


deu 

~dt 


exx  ■  Vfi{x^  )-ex2  ■  V f2{x^  )  +  ex^  ■  V h{x^  ) 


(A.2) 


1  u^(t  X ) 

Defining  the  eorresponding  nominal  angle  0^(t,  x  )  =  tan~^(  and  re-arranging  the  terms 

'^X  ) 

of  Eq.  (A.3),  we  get  the  linearized  version  of  the  geometrie  heat  equation  around  a  nominal  value: 


^  «  AGmr,Mt,x) 

Ot 

=  sin^(0”(f,  X  ))  Uxx{t,  X  )  —  sin(20”(f,  x  ))  Uxy{t,  x  )  +  cos^(6”(f,  x  ))  Uyy{t,  x  ) 
+  c  (-Uy(t,x)ux{t,x)+Ux{t,x)uy(t,x)),  (A.4) 


where  c  =  [sin(20”)«^  -  u^y) 


cos{20‘^)2Uxy  ■ 
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Appendix  B 


Derivation  of  DDEs  for  vertex  motion 


Let  us  define  C  y  for  vertex  Vi  by 


Cy{p,V) 


{p- 

-  (*  -  l))e 

for 

i  —  1  <  p  <  i 

(1- 

-  {P  -  *))e 

for 

i  <  p  <  i  +  1 

0 

for 

\p  —  i\>l 

where  e  =  {ex,  e  y)^  denote  the  standard  bases  for  .  We  now  write  our  energy  as  a  funetion  of 
the  vertiees  V  :  E{V)  =  f^(E,  N)  HCpUdp  =  J^iF,  JCp)  dp,  and  eompute  its  partial  derivative 
with  respeet  to  one  of  the  vertex  eoordinates  v  where  either  v  =  Xi  or  v  =  yi  for  some  1  <  i  <  n, 
as  Ey{V)  =  Jq  f{Cy,  JCp)  dp.  Substituting  the  partieular  forms  of  C ,  C  p  (Eq.(4.4)),  and  C  y 
(the  latter  denoting  either  C  xi  or  C  y.  for  some  vertex  i)  into  this  expression  yields 


Ev=  j^f{L{p-\j)\,V  [pj  5  [pj+l))  V,  J{^  \p\-\-l 

n-1  „i 

=  {Cv{p  +  k),J{Vk+i-Vk))dp 

=  Jj{L{p,Vi.i,Vi)){Cy{p+i-l),JiVi-Vi.i))  dp 
+  IjiLiP,  Vi,  Vi+i)){Cyip  +  i),  JiVi+i  -Vi))dp 
=  (e ,  J{V  i  pf  {Lip,  Vi.i,V  i))  dp 

+  {e  ,  J{Vi+i  -Vi))  l\l  -  p)f  {L{p,  Vi,V  i+i))dp 
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where  e  =  e a;  if  u  =  rr j  or  e  =  Sy  v  =  yi.  If  we  introduee  a  time  variable  t  and  evolve 
both  coordinates  Xi  and  y*  in  the  gradient  directions  given  above,  and  denoting  the  corresponding 


J{V  i  —  V  j_i)  =  N  we  obtain  the  following  gradient  flow  for  the  vertex  Vi 


dVj 

dt 


Vi.i,Vi))dpNi,i_i+  [  (l-p)f{L(p,Vi,Vi+i))dpN 

Jo 
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Appendix  C 


Derivation  of  the  Gradient  Flow  for 
Varying  Priors  proportional  to  areas 


Rewriting  here  Eq.(4.19),  (note  that  A  =  \Ru\  +  |-R«|), 


dC 


=  VJSa,2  = 


E  (2(«y  -  v,)(Vui  -Vvj))N, 


i=i 

1  \Rv 


l\R^ 


2  A2 


2  A2 


(C.l) 


it  may  be  rewritten  as  follows: 


^  f^(u  V  )  f^l^uHRvl(Gj(I)-uj)  ^  2IRYIRY(Gj(I)-vj 


dt  2A2 


V 


i=i 

+  (|^«|  -  \Ru\)iUj  -  Vj] 


\Rn 


\Rv 


2A2 

2A2 

2A2 


'^{Uj  -  Vj)  {2\Ry\{Gj{I)  -  Uj)  +  2\Ru\iGj{I)  -  vj)  +  (|i?„|  -  \Ru\)iuj  -  vj)) 


i=i 


^^{uj  —  Vj)  {2\Ry\Gj{I)  —  2\Ry\uj  +  2\Ru\Gj{I)  —  2\Ru\vj  +  (|-R?)|  —  \Ru\){uj  —  Vj)) 


i=i 


Vj)  {2{\RY\  ~\~  \Rv\)G ji^I) 


i=i 


^  -  Vj)iGjiI)  -  Uj  +  Gj  -  Vj) 


i=i 


(C.2) 


i=i 
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It  can  be  observed  that  Eq.(C.2)  is  the  gradient  descent  flow  for  the  following  energy  functional 


dxdy  +  /  (Gj(I)  —  VjYdxdy, 


(C.3) 


which  is  a  generalized  form  of  the  energy  functional  proposed  by  Chan  and  Vese  [29]. 
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Appendix  D 


Derivation  of  the  Electric  Field  at  a 
point  exerted  by  a  line  charge 


The  integral  we  want  to  evaluate  is  given  by 

rb 


Eabix')  =  f 

Ja 


X  —  X 


rXdx 


\x'  —  x\ 


=  kr  L 


Jo  ll(^' 


x'  —  a)  +  t{a  —  b  ] 


—  a)  +  t{a  —  b  ] 


:dt. 


(D.l) 

(D.2) 


By  making  ehange  of  variables,  u  =  x'  —  a,  and  v  =  a  —  b,  the  integral  above  ean  be  written  as 

Eab=  C  (D.3) 

Jo  \\u+tv\\^  Jo 

where  we  denoted  the  integrand  by  a  veetor  f'{t),  a  funetion  of  t,  whieh  is  the  derivative  of  the 
solution  /  (t)  w.r.t  t.  We  expeet  the  solution  /  (t),  to  this  integral  equation,  to  be  of  the  form 


fit) 


{au  +  bv)  +  t{cu  +  dv  ) 
lilt  +  II 


(D.4) 


where  a,  b,  c,  and  d  are  sealars  to  be  estimated.  Taking  the  derivative  of  /  (t)  w.r.t.  t,  we  obtain 
(note  that  in  the  following,  for  simplieity  we  denote  the  inner  produet  of  a  veetor  by  itself  v  by 
V  and  the  inner  produet  by  •) 


fit) 


{cu  +  dv){v‘^  +  2tu  V  +  t^v  —  [(ait  +  bv)  +  t{cu  +  dv)]{uv  +  fu 

I  |it  +  fu  I P 


(D.5) 
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The  numerator  above  ean  henee  be  equated  to  the  numerator  in  the  original  integral,  i.e.  u  +  tv  . 
This  leads  to  the  following  two  equalities: 

u  =  u{cu‘^  —  a{uv  ))  +  V  {du"^  —  b{uv  )) 

V  =  u{c{uv)  —  av"^)  +  V  {d{uv)  —  bv"^)  (D-6) 


Equating  the  eorresponding  eoeffieients  of  u  and  w  to  0  and  1  as  required  leads  to  the  following 
system  of  linear  equations  in  four  unknowns: 


—  u  V 

0 

u  ^ 

0 

a 

1 

0 

—  u  V 

0 

b 

0 

—  V  ^ 

0 

U  V 

0 

c 

0 

0 

—  V  ^ 

0 

U  V 

d 

1 

Carrying  out  Gaussian  elimination  on  this  system  of  equations  leads  to  a  unique  solution 


a 

—  u  V 

b 

_  1 

~  {uv^-u’^V^ 

9 

c 

—  V  ^ 

d 

U  V 

Substituting  the  sealar  eoeffieients  a,  b,  c,  and  d  into  the  Eq.(D.4),  the  solution  to  the  integral  ean  be 


written  as 


fit) 


—  {uv)u  +u‘^v  +t{—v‘^u  +{uv)v) 

{{uv)‘^  —  u‘^v‘^)\\u  +  II 


1 

0 


(D.7) 


The  resulting  eleetrie  field  af  x exerted  by  fhe  line  eharge  exlending  from  a  lo  6 ,  is  Ihen  given  by 


kcL 

1  {uv  +  v‘^)u  —  +  U  V  )v 

{uv)u  —  u^v  1 

{u‘^v‘^  —  {uv  )2) 

\  tt  +  w 

ll^ll  i 

(D.8) 


Noting  lhal  (see  Eigure  4.8) 


/ 

def 

u 

=  X 

-  a  =  Xa 

U  +  V 

/ 

,  def 

=  X 

-b  =  Xb 

,  def 

V 

=  a 

-  b  =  X  b  —  X 
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and  substituting  the  above  veetors  into  Eq.(D.8),  we  obtain 


{Xg-  {xb-Xa)  +  \\xb-  Xa\\‘^)Xa  -  {\\x  a\\‘^  +  X  g  ■  {x  b  -  X  a)){x  b  -  X  g) 

W^bW 


X  g  ■  {x  b  -  X  g)x 

a  1 

\^a\ 

P(a:6  -  aJa)  1 

[  kcL 

1 

\Xa\ 

1  J 

f  1 

\Xg\\‘^\\Xb  -Xg\ 

1“^  -  {Xg-  {Xb  -Xg)f 

(D.9) 


whieh  ean  be  simplified  to  Eq.(4.31)  rewritten  here  for  eonvenienee 
Egb(x')  = 


he  L 

1 

'\\Xb\\^Xg 

-  {Xg 

■Xb)Xb  , 

\\Xg\\^Xb  -  (x  g  ■  X  b)x  g 

||a:ain|a36|p  -  {Xg 

■xbY  \ 

\xb\\ 

+  ■ 

Hawaii 

(D.IO) 
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